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Chapter One 


Introduction and Summary 


1.1 Background 

On July 20, 1989, the 20^ anniversary of the first Apollo moon landing. 
President George Bush challenged the nation to undertake an ambitious 
course of human space exploration. After establishing a manned presence in 
earth orbit with the Space Station Freedom in the 1990's, the President 
proposed that the U.S. return to the moon, and return to stay. From this 
lunar basing point, the U.S. could continue human exploration of our solar 
system by undertaking a manned mission to Mars. 

The establishment of a lunar base will result in significant lunar traffic 
to build, supply, and resupply this facility. This increased traffic will require a 
lunar navigation system. As the nation prepares its return to the moon, it 
will have to decide whether this navigation system should be earth-based, 
vehicle-based or lunar-based. 

Our initial voyages to the moon primarily depended upon earth-based 
navigation systems, although the manned missions had some on-board 
capability. An earth-based method could be adopted for future lunar travel, 
but NASA's Deep Space Network (DSN) tracking is manpower intensive, 
costly, and is not suitable for high traffic rates. Additionally, earth-based 
navigation can only track vehicles on the lunar near side. Adopting an earth- 
based navigation system would be an impractical stepping stone for human 
exploration of Mars and the solar system. 

A vehicle-based navigation system could be developed to support 
future lunar traffic. Vehicle-based navigation has become practical because of 
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advances in inertial navigation equipment and on-board computing 
capabilities. Inertial navigation would be limited by our knowledge of the cis- 
lunar environment, principally our knowledge of the moon's gravity field. 

The limitations of earth-based navigation and the high accuracy 
requirements of certain mission phases (principally landing) may require the 
development of a lunar-based navigation system. A simple lunar navigation 
system similar to the earth's Global Positioning System (GPS) would handle 
high traffic rates and would provide accurate lunar far-side navigation. If 
high lunar traffic rates are achieved, then a lunar-based system could provide 
a higher accuracy system alternative to vehicle-based navigation. 

Many low-altitude mission phases will require accurate knowledge of 
the moon's gravitational field, especially its far-side characteristics. Spherical 
harmonic models are typically used to model the gravitational field near 
celestial bodies. Finite expansion spherical harmonic models, however, do 
not accurately model the low altitude gravity field of the moon. This is 
because the moon's gravitational field contains significant anomalies, 
discovered in the late 1960's by scientists at NASA's Jet Propulsion Laboratory 
(JPL) [35], making such models inefficient. From earth-based lunar tracking 
data, the scientists developed a lunar gravity model. This model was then 
compared to topographical images of the moon and revealed mass 
concentrations around the ringed maria. These mass concentrations or 
"mascons" exhibit very high frequency gravitational behavior and therefore 
require a very high number of terms in the spherical harmonic expansion to 
model this behavior. Since their discovery, scientists have postulated 
different models to account for the lunar mascon phenomenon, since 
expanding the spherical harmonic model to high degree and order was 
computationally impractical. Chapter Two surveys the spherical harmonic 
and several other gravitational field modeling techniques in more detail. 


1.2 Motivation 

Since the real lunar gravitational field is difficult to accurately model, 
this thesis will study the implications of modeling errors. An inaccurate 
model of the lunar gravity field will result in the growth of navigation errors. 
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Mismodeled acceleration forces result in both velocity and position errors. 
Since gravitational acceleration is a function of position, position errors will 
lead to increased acceleration errors, further increasing the velocity and 
position errors. This error propagation may or may not be critical depending 
upon the magnitude of the errors, the navigation system's ability to measure 
them, and the mission phase accuracy requirements. 

Specifically, an inaccurate gravity model will significantly affect any 
landing maneuvers with strict accuracy requirements. Unmanned cargo 
missions to resupply a lunar base will be particularly vulnerable to errors 
from a mismodeled gravity field. Since there is no appreciable lunar 
atmosphere, gravity forces dominate a vehicle's descent to the moon's 
surface. Since the force of gravity is inversely proportional to the square of 
distance, navigation errors due to a mismodeled gravity field increase as the 
vehicle descends to the surface. Avoiding unacceptable landing errors will 
depend upon an accurate determination of the lunar gravity field. 

The scientific community is also interested in developing a more 
precise model of the lunar gravity field. A better model can improve 
knowledge of the moon's composition and internal structure. Models of 
different elements, their densities, and their distribution within the moon's 
interior could be developed to match the observed gravitational field. 
Gravitational models may also help to determine the selenological thermal 
and tectonic history. The discovery of lunar mascons has also led to scientific 
speculation about how mass concentrations formed in these shallow seas. 
The scientific community hopes that a better understanding of the 
gravitational field around the ringed maria and other lunar surface features 
will help to determine the origin of these features [2, 35]. 

The purpose of this thesis is to determine the feasibility of using a 
spherical harmonic lunar gravitational model, based on observations of a 
near-circular polar satellite, to predict low altitude lunar orbits globally. 
Rather than attempting to develop a more precise lunar gravitational field 
model, this thesis investigates measurement types and satellite orbits that can 
be used to develop gravity field models. Each measurement type will have 
advantages and disadvantages in terms of cost, schedule, and accuracy. This 
thesis investigates each different method's ability to estimate a lunar 
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gravitational potential model and the model's accuracy in predicting future 
lunar orbits. By analyzing the capabilities of different sensing methods, this 
investigation will allow NASA to plan unmanned lunar precursor missions 
to extract the best lunar gravity field information. 


1.3 Initial Lunar Gravitational Sensing Method 

Current lunar gravitation field models are based upon earth-based 
tracking data from the Lunar Orbiter program of the 1960’s (Figure 1.3-1). 
These unmanned Apollo precursor missions provided photographic imaging 
and gravitational mapping of the moon. Apollo lunar navigation was based 
upon Lunar Orbiter's gravity field mapping. In addition to the Lunar Orbiter 
missions, tracking data from Apollo missions and some Soviet lunar 
missions are included in current gravity field models [3,12,19, 33, 41, 47]. 

Lunar Orbiter gravitational mapping missions utilized Doppler 
measurements of radio tracking signals. Lunar Orbiter spacecraft were tracked 
by NASA’s Deep Space Network (DSN) across the near side of the moon [36]. 
A DSN tracking station sent a continuous wave S-Band frequency to the 
spacecraft. The spacecraft received this Doppler-shifted signal and re¬ 
transmitted it to earth. The tracking station received this signal, Doppler- 
shifted once again in frequency. The tracking station used this signal to 
calculate the relative velocity between the spacecraft and tracking station. The 
relative velocity observed was then combined with position tracking data to 
estimate the lunar gravity field using methods similar to those discussed in 
Chapter Five. 



Figure 1.3-1: Lunar Orbiter Earth-Moon Geometry 


Current lunar gravitation field model accuracy is limited by the 
amount of lunar orbital tracking data available. The Lunar Orbiter and 
Apollo missions were mostly low inclination missions. Since most missions 
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flew about the lunar equator, derived gravity field models emphasize the 
effects of anomalies near the equator. Only a fraction of the moon's surface 
was covered by these missions [32], and therefore gravity field data is lacking 
for lunar polar regions. Additionally, this earth-based gravity mapping 
method was limited to observations of near-side lunar spacecraft passes. 
Gravitational disturbances on the far side were only determined by their 
integrated effects on satellite position and velocity from the end of one near¬ 
side pass to the beginning of the next pass. Thus current lunar gravitational 
field models do not provide very meaningful information about the lunar far 
side. 


1.4 Proposed Lunar Gravitational Sensing Methods 

Any future lunar gravitational field sensing system will have to greatly 
improve our knowledge of the moon's gravitational field to justify the 
mission's cost. To achieve this improvement in accuracy, the system will 
have to address the current model's limitations. The motion of an orbiting 
body should be sensed without any orbital maneuvering which disturbs the 
estimation solution. 1 Thus it is desirable to select orbits which are stable for at 
least one lunar orbit to avoid re-boost maneuvering. A high inclination, 
preferably polar, lunar orbiter would allow observations of satellite 
accelerations over the moon's entire surface as the moon rotates under the 
orbital plane. A dual orbiter sensing scheme would be better because it would 
allow lunar far-side accelerations to be observed. Better still would be a 
sensing scheme observing the motion of several satellites in different 
inclinations than those available in Apollo-era lunar missions. 

NASA is considering two different sensing schemes. NASA's Jet 
Propulsion Laboratory has proposed a dual orbiter scheme which uses radio- 
based Doppler observations to sense the moon's gravitational field effects [40]. 
NASA's Goddard Space Flight Center has proposed a co-orbital scheme which 


1 Gravitational accelerations experienced by an orbiter are not measured directly. Methods to 
measure the accelerations due to gravity therefore use external observations of a body's motion. 
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uses a laser-based system which makes both ranging and Doppler 
observations [2], 

The dual orbiter sensing scheme uses a low altitude, circular polar 
satellite and a high apolune polar elliptical satellite (Figure 1.4-1). The 
elliptical satellite orbit is positioned so that apolune initially occurs on the far 
side of the moon. This increases the duration of lunar far-side viewing. The 
orbit is also skewed such that apolune is outside of the earth occultation zone 
which increases the time that the elliptical satellite is within the line of sight 
of earth tracking stations. From these orbits the relative velocity between the 
two spacecraft can be measured using either the bent pipe or the satellite 
bounce methods described below. 



The bent pipe method uses a four-way coherent Doppler scheme in 
which a high frequency is generated by an atomic clock at a DSN tracking site 
and transmitted to the elliptical "viewing" satellite. The "viewing" satellite 
uses the Doppler-shifted received signal to generate a lower frequency signal 
which it transmits to the circular "gravity sensing" spacecraft. This spacecraft 
receives the Doppler-shifted signal and re-transmits it to the "viewing" 
satellite. The "viewing" satellite modulates the received Doppler-shifted 
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frequency signal from the "sensing" spacecraft onto the frequency signal 
received from the tracking station and transmits this signal back to earth. The 
signal received by the tracking station is then processed to retrieve both the 
relative velocities between the tracking station and "viewing" satellite and 
between the "viewing" and "gravity sensing" satellites. 

The satellite bounce method uses a two-way coherent Doppler scheme 
between the two spacecraft. The circular "gravity sensing" satellite generates a 
continuous wave frequency signal for transmission to the elliptical "viewing" 
satellite. This satellite shifts the signal's frequency for transmission back to 
the first one. The receiving spacecraft extracts the Doppler shift from the 
signal, records it and transmits it to earth when in view of an earth tracking 
station. Coherent Doppler links between earth tracking stations and either 
satellite are used to aid in the estimation of the lunar gravitational field. 



NASA's Goddard Space Flight Center has proposed a co-orbital scheme 
in which a satellite in a circular polar orbit ejects a subsatellite in the same 
orbit (Figure 1.4-2). Both spacecraft are affected by lunar gravitational 
perturbations, so both are "sensing" vehicles and a laser system measures 
their relative motion. The satellite contains the sensing equipment and the 
subsatellite is a passive reflector for the satellite's emitted laser beams. The 
satellite's laser transmits a light beam toward the subsatellite. This beam 
reflects off of the subsatellite back to the satellite. A high accuracy ranging 
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measurement between the satellite and subsatellite is made by determining 
the travel time of a transmitted sub-carrier pulse signal- Lesser accuracy 
Doppler relative velocity measurements are made from the frequency 
shifting of the transmitted laser signal. Observation data is stored and 
transmitted to earth at regular intervals. 


1.5 Simulation Tools 

The primary tool used to accomplish the goals of this thesis is the 
Planetary Ephemeris Program (PEP), a FORTRAN computer program 
obtained from the Smithsonian Astrophysical Observatory (SAO) and 
executed on Sun workstations at the Charles Stark Draper Laboratory. The 
Smithsonian version of PEP has most of the capabilities needed for the 
analyses of this thesis. Modifications, coded at Draper Laboratory, have 
augmented its capabilities for this thesis research. In addition, auxiliary 
software has been developed to analyze lunar orbits, the observation schemes, 
navigation uncertainties, and estimated lunar gravitational fields. 

Given the description of a body's gravitational field, PEP can 
numerically integrate a satellite's motion about that body. For this thesis, PEP 
was modified to accommodate a point mass (mascon) gravity model in 
addition to the spherical harmonic model. The techniques used in PEP for 
numerically integrating the differential equations of motion for a lunar 
satellite are described in Chapter Three. This chapter also describes the 
methods PEP uses to calculate the partial derivatives of the motion with 
respect to orbital initial conditions, gravity harmonic coefficients, and other 
parameters. The mascon gravity model modifications coded in PEP are also 
covered in Chapter Three. 

The Planetary Ephemeris Program (PEP) can also process the 
astronomical observations generated by the various lunar gravitational 
sensing methods. The many different types of observations that can be 
processed in PEP are described in Chapter Four. This chapter also discusses 
how PEP generates simulated observations with a truth model (mascon and 
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spherical harmonic expansion) and then estimates theoretical values of the 
observations for another model (spherical harmonic expansion) to fit those 
"truth model" observations. Least squares maximum likelihood estimation 
and prediction uncertainty propagation techniques are described in Chapter 
Five. For this thesis, noise has not been included in the satellite dynamics. 
Kalman filter and system identification techniques may be needed when 
processing real observations, because of the noise due to radiation pressure, 
gas leakage, and other unmodeled forces. These estimation techniques are 
also described in Chapter Five. 


1.6 Methodology 

The focus of this thesis is the estimation of a lunar gravity field model 
based on various measurement types and/or orbital geometries. Each unique 
measurement type and orbital geometry combination will be refered to as a 
sensing scheme. The standard earth-based orbiter state sensing scheme and 
the proposed dual orbiter bent pipe scheme are analyzed in-depth. 
Additionally, the co-orbital laser ranging scheme, a non-coplanar bent pipe 
scheme, and an earth-based interferometric observation scheme are 
investigated to determine whether any of these schemes can reduce the 
parameter correlation's observed during gravitational parameter estimation. 

Since the true lunar gravitational field is not precisely known, a 
"truth" model was developed and used for this investigation, as discussed in 
Section 6.4. This "truth" model combines Bills and Ferrari’s 16 X 16 lunar 
harmonic model [12] up to degree and order five, along with 78 point masses 
distributed below the lunar surface to simulate the behavior of mass 
concentrations. This truth model was used to simulate observations for the 
various sensing schemes. 

For each sensing scheme, the coefficients in spherical harmonic fit 
models of degree and order 8 and 12 were estimated to optimally represent 
the "true" gravity field model by fitting to the "truth" model observations. 
Using first guesses for the fit model coefficients and satellite initial osculating 
orbital elements, the equations of motion and the equations for the partial 
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derivatives of the motion with respect to these quantities were numerically 
integrated. From these numerical integrations, theoretical values of the 
observations were determined. The observation residuals (difference 
between the "truth" model observations and fit model's theoretical 
observation values) and the observation partial derivatives were computed 
and used to obtain parameter adjustments to the first guesses for the fit model 
coefficients and satellite initial osculating orbital elements. This process was 
then repeated until either the method converged upon a solution or no 
solution could be determined. Chapter Six discusses the implementation of 
the estimation process for various sensing methods and certain test cases. 

The estimated lunar gravitational field models were analyzed to 
determine their accuracy relative to the "truth". Navigation errors 
propagated in one lunar revolution were used to determine the estimated 
gravitational field model’s accuracy. Two types of lunar orbits were analyzed: 
a 15° inclination, 100 km altitude, near-circular orbit and a lunar landing 
from a 5° inclination, 200 km near-circular parking orbit. These orbits were 
propagated using both the "truth" and estimated models. The analysis of the 
estimated gravitational field models is discussed in detail in Chapter Seven. 
In addition, this chapter discusses the attempts made to break the high 
parameter correlations which were discovered during the estimation process. 

For the near-circular orbit, the position and velocity errors between the 
"truth" and estimated gravity fields were used to quantify the estimated 
model’s accuracy and thus the sensing method's capability. The gravitational 
parameter covariance matrix, determined during the estimation process, was 
used with the estimated field's orbit propagation with partial derivatives to 
predict position and velocity uncertainties. These predicted uncertainties 
were then compared to the true state errors between the two models' orbits. 
In real gravity field missions, the true gravity field will not be available for 
comparison with the estimate, so it is useful to understand the relation 
between these two analyses. 

For the lunar landing maneuver, the estimated gravity field model was 
used to determine the deorbit burn and the selenographic position for 
Powered Descent Initiation (PDI). The spacecraft's circular parking orbit was 
numerically integrated until the appropriate time for the deorbit burn. The 
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spacecraft's velocity at burn time was then adjusted to simulate the deorbit 
maneuver and the spacecraft's elliptical transfer orbit was then integrated. 
The spacecraft's position upon reaching the powered descent stage of the 
mission was then determined. The "true" PDI point was then compared to 
the estimated model target PDI point to determine the model's accuracy for 
planning lunar landing maneuvers. 


1.7 Summary of Results 

The analyses of the 8x8 estimated fit models clearly demonstrated that 
lunar far-side observations are required in the accurate estimation of the 
lunar gravity field. For both the lunar landing maneuver and the satellite 
state uncertainty prediction, the observation technique which included lunar 
far-side observations in addition to earth-based near-side observations 
produced a much more accurate lunar gravity fit model. This estimated 
model planned a lunar deorbit maneuver 4.3 times more accurately than the 
model based on earth-based observations alone. The earth-based observation 
fit model also predicted state uncertainties four times larger than its 
counterpart. The earth- and satellite-based fit model produced state errors for 
the single orbit that were again roughly a quarter of the earth-based fit 
model's errors. 

The lunar navigation analyses also demonstrate that the eighth degree 
and order spherical harmonic expansion fit models were unable to adequately 
model the lunar mascons included in the lunar gravitational "truth" model. 
In the best case, the 8 x 8 fit model predicted single orbit ahead uncertainties of 
close to three quarters of a kilometer in position and one meter per second in 
velocity. The orbit's actual state errors were closer to three kilometers in 
position and two and a half meters per second in velocity. Additionally, the 
best fit model produced a fifty-six kilometer position error for the lunar 
deorbit maneuver. These results are discussed in further detail in Chapter 
Eight, which also recommends several subjects related to this thesis which 
deserve further study. 
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Chapter Two 


Gravity Field Models 


2.1 LaPlace's Equation and the External Gravitational Field 


According to Newton's law of gravitation, two particles attract each 
other with a force, acting along the line joining them, which is proportional 
to the product of their masses and inversely proportional to the square of the 
distance between them [10, p. 95]. From this law, the gravitational force of a 
body of mass m 2 acting upon a body of mass mi can be mathematically 
represented by the formula: 




X 2 ~ T \ 


( 2 . 1 - 1 ) 


where G is the universal gravitational constant, and r, and r 2 are the position 
vectors of bodies one and two respectively. 


Unfortunately, this formula is only appropriate if the two bodies are 
point masses, or behave as them. Such is the case for spherical bodies if 
density is a function of radius from the center only. The point mass model 
also provides an accurate representation of the gravitational attraction for 
widely separated bodies. In the limit of large distances, gravitational bodies 
tend to look like point masses so that mass distribution becomes 
unimportant. 


For many practical applications, the attracting body cannot be modeled 
as a point mass and a different mathematical model must be developed. For 
the case in which the attracted body is small compared to the attracting body 
and the attracting body is an arbitrary distribution of mass with a finite 
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dimension, the force on a mass m located at position vector r in Figure 2.1-1 
produced by an element of mass dM with position vector R is 

(f-K) 

dF(r) = -Gm dM - zj 

f-R 



Integration over the entire volume of the distributed body will produce the 
gravitational force on the mass m of the attracting body of total mass M. This 
force can then be represented by a scalar potential U, such that the 
gravitational force on a body located outside of the attracting body may be 
obtained as the gradient of the scalar potential, or 

F(r)=-mVU(f). (2.1-3) 

There is a sign convention discrepancy between some of the references used 
and PEP documentation [8]. According to Kaula [24], physicists define the 
gradient of the potential field as in Equation (2.1-3), whereas astronomers 
define the same gradient with a change in sign. The formulas in this thesis 
follow the former convention to agree with PEP documentation and software 
coding modifications explained herein. 


The scalar gravitational potential U can be written in terms of the rectangular 
coordinate system of Figure 2.1-1 and the mass density p as [14] 






[(*-£) 2 +( y _f )) 2 +( z - f ) 2 


K 


d&Tidt 


(2.1-4) 
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The Laplacian of the scalar potential U in rectangular coordinates is 


V 2 U = 


<? 2 U d 2 U d 2 U 


(2.1-5) 


dx 2 dy 2 dz 2 

Taking the required partial derivatives of the potential U from (2.1-4) yields 

-3 


V 2 U = -GM 


(*-'5) 2 +(y-'i) 2 +( z -0 

(x-4f +(y-v) 1 +(z-Q 


i Yi 


+ 


+(y-i ) 2 +(z-Cf 


(2.1-6) 


Since the bracketed terms cancel, this reduces to 

V 2 U = 0 


(2.1-7) 


This relationship is known as Laplace's equation and applies at all points 
outside of the distributed attracting mass. Its solutions are called harmonic 
functions. Any scalar function, U, which satisfies Laplace's equation and the 
far-field boundary condition that the potential approaches 0 as 1 /r can be used 
to descibe the gravitational field about some distributed mass. If U is defined 
with sufficient flexibility, i.e. an infinite number of orthogonal terms with 
undefined constant coefficients, then U can be tailored to describe the 
gravitational field outside of any arbitrarily distributed mass. The above 
method of deriving (2.1.7) is based upon the method used by Kaula [24], Battin 
[10], and Comfort [14]. 


2.2 Spherical Harmonic Expansion for the Gravitational Potential 

The most common gravitational potential model is the spherical 
harmonic expansion. This expansion can be derived by solving Laplace's 
equation in spherical harmonic coordinates. First, the rectangular coordinate 
system of Figure 2.1-1 is converted to spherical coordinates through the 
transformation below, which is depicted in Figure 2.2-1. If the center of the 
distributed mass is selected as the origin of this coordinate system, the 
expansion is simplified. 
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X = rcos0cos0 
y = rsin0cos0 
z = rsin 0 


0 < r < °o 

0 < 9 < 2 k 
-%<<!><% 



( 2 . 2 - 1 ) 


Figure 2.2-1: Spherical Coordinate System 
In spherical coordinates, Laplace's equation becomes [24] 


r 2 V 2 U = f- 

dr 


d( 2 d U V 1 


V 


dr 


cos 0 t?0 


d ( .dU 
cos0— 
d<t> J 


1 <9 2 U 

cos 2 0 r?0 2 


s 0 (2.2-2) 


One common method of deriving the solution to U in spherical 
coordinates is through the method of separation of variables, maintaining the 
boundary condition on r. This method can be found in Kaula [24] or Comfort 
[14] and leads to the solution 

U(r / 0,0) = ^-^ r ^P nm (sin0)[C Mm cosm0 + S wrn sinm0] (2.2-3) 

«=0 r m =0 

where C nm and S nm are constant coefficients and P n m are the generalized 
Legendre functions of degree n and order m. Equation (2.2-3) is the complete 
real solution of Laplace's equation in spherical coordinates. There are other 
solutions which are physically inadmissible since they cause the potential to 
become infinite at 0 = k/ 2, 3n/2. These solutions involve the associated 
Legendre functions of the second kind. Another method of deriving the 
solution to U is to expand the denominator of (2.1-4) in Legendre functions. 
This is the procedure carried out by Battin [10] and Croopnick [18]. 

Equation (2.2-3) is the generalized solution of Laplace's equation and 
must be modified for the gravitational case. Laplace's equation holds true 
everywhere outside the surface of the gravitational body. The harmonic 
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expansion, however, is convergent only from pj to <*> in the radial direction, 
where p/ is the body's (lunar) maximum radius, usually the equatorial radius. 
Furthermore, pi can be used to non-dimensionalize the C nm and S nm 
coefficients as shown in Equation (2.2-4). Secondly, U is multiplied by the 
gravitational constant and the total mass of the attracting body to satisfy 
Equation (2.1-3). The n=0, m=0 term is the point mass approximation for the 
attracting body and can be separated from the summation for ease of 
computation. The S nm coefficient is meaningless for m=0 because sin(m0) is 
always zero. Separating the m=0 terms simplifies the summation. The 
P n0 (sin<t>) Legendre functions then become the standard P n (sin<])) Legendre 
polynomials. When the m=0 terms are separated, the C n o terms are 
commonly replaced by J n terms, where J n = -C n o- After the previous 
modifications, the spherical harmonic gravitational potential, U, becomes 


u M,0) = 


-GM 


i-yjJ —1 P„(sin 0 ) + 
Zt vr; 


, n n 


X — ^Pnrn^^nm^sme + S^smmG] 
n =1 '• ^ ' m =1 

(2.2-4) 


Additionally, if the attracting body's center of mass is the origin of the 
coordinate system, then the first degree terms are identically zero [22], 

Ji=0, C„=0, Sii = 0 (2.2-5) 

and the summation limits can begin with the second degree terms. 


When dealing with large degree models, it is common to normalize 
the Legendre functions, and therefore adjust the tesseral coefficients as well. 
The Legendre functions are normalized to satisfy the equation 


Jj* (jP nm (sin <p) cosm #) 2 dOd(f> = 4n 
JJ (P„ m (sm 0 )smm 0 ) 2 d 0 d 0 = 4 n 


0<6 <2k 

\-y 2 <(t)<y 2 


(2.2-6) 


This leads to the following relationship between the unnormalized and 
normalized Legendre functions and the inverse relationship between the 
unnormalized and normalized coefficients: 


35 


LUNAR GRAVITATIONAL FIELD ESTIMATION AND SATELLITE ORBIT PREDICTION 


Pn= V2^IP„ 
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2{2n + l)(n - m)l 
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(m = 0 ) 
(m> 0) 


(m = 0 ) 
(m> 0) 


(2.2-7) 


( 2 . 2 - 8 ) 


The Planetary Ephemeris Program (PEP) was modified to allow 
normalized Legendre functions to be used for higher degree harmonic 
models. Normalized coefficients are better conditioned for PEP's floating 
point parameter estimation algorithms, especially for higher degree models. 
Appendix A discusses the Legendre polynomials, the generalized Legendre 
functions and the recursive formulas developed for their implementation. 
Although use of the normalized Legendre functions was considered necessary 
for this thesis, this scaling was not considered necessary for the Legendre 
polynomials. Because of this inconsistency in scaling, the equations coded in 
the software use normalized tesseral coefficients and unnormalized zonal 
coefficients. 


After accounting for the above-mentioned modifications, the spherical 
harmonic potential equation (2.2-4) becomes 


U(r,<M) = 


-GM 


p »( sin 0) 

n=2 Y ' 


+ 


ln=2 


. ft n 


XlT") J,[C nm cosme + S nm smrne]P nm (sm<l>) 

(2.2-9) 


m =1 


The Planetary Ephemeris Program uses this equation with finite summation 
upper limits. The coding in PEP separates the central body term from the 
harmonic summations. 
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2.2.1 Relation between Spherical Harmonics and Moments of Inertia 


The moments of inertia about a solid body are closely related to the 
gravitational potential for a distributed mass. An expansion of Equation (2.1- 
4) for the gravitational potential contains the terms for the moments of 
inertia. 


«=If/o» 2 < 

| dm 

kv=!H^ dm 

body 


body 

„=J 

) dm 


body 


body 

« = jjj ^ 2 + 7 ^ 

body 

) dm 

I,{ = JjJ VC d™ 

body 


Equations (2.2.1-1) are valid for any arbitrary rectangular coordinate 
system (£,t|,£) originating at the body's center of mass. The spherical 
harmonic potential function U was developed from the (r,<)>,0) coordinate 
system of Figure 2.2-1 which corresponds to the (x,y,z) rectangular coordinate 
system of Figure 2.1-1. Converting the moments of inertia to this (x,y,z) 
system, the second degree coefficients in the gravitational potential function 
satisfy the following relations [24, 30, 32,44]. 


^20 “ 
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> 22 ' 




2 pfM 


Since the lunar moments of inertia can be obtained by observing the 
physical librations of the moon, the second degree coefficients can be 
determined without sending spacecraft to the moon. If the (x,y,z) axes are 
principle axes, the products of inertia will be zero, resulting in 

C 2 i = 0, S 2 i = 0, S 22 = 0 (2.2.1-3) 


In the lunar case, the principle axes coincide with the x axis pointing 
towards the earth and the z axis pointing along the axis of rotation. 
Unfortunately, equations (2.2.1-3) do not hold exactly with the inertial 
reference frame used by PEP in this thesis. Since the lunar moments of 
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inertia can be experimentally determined from earth, it is possible to 
simultaneously process lunar rotation observations with satellite 
observations such as those simulated in this theses. In this manner, 
equations (2.2.1-3) and the non-zero equations in (2.2.1-2) could be assumed to 
hold exactly. Processing lunar laser corner reflector observations [13, 17] 
would provide estimates of the second degree coefficients with a very high 
degree of confidence and the lunar satellite observations could be 
simultaneously processed for estimates of the higher degree coefficients. 


2.2.2 Limitations of the Spherical Harmonic Expansion 

Although the spherical harmonic expansion is frequently used to 
describe the gravitational potential of a distributed mass, it is not an efficient 
model for all uses. Since the model is based upon spherical coordinates, it 
produces an uneven resolution of coverage from the equator to the poles. 
This uneven resolution can lead to modeling inefficiencies. 



Figure 2.2.2-1: Zonal, Tesseral, and Sectorial Harmonic Patterns 


The spherical harmonic model breaks the sphere up into zonal, 
tesseral, and sectorial patches (Figure 2.2.2-1). The patches are separated by 
lines where the terms are identically zero. On one side of the line, the term 
will be positive and on the other it will be negative. The zonal, J n/ terms, are 
independent of longitude, 0, and dissect the globe of interest into n+1 bands 
along lines of latitude. The tesseral terms, C nm and S nm , dissect the globe of 
interest into patches of both latitude and longitude. A tesseral term will have 
n-m+1 sections of latitude and 2m sections along lines of longitude. The 
sectorial terms, C nn and S nn , are independent of latitude, <|>, and dissect the 
globe of interest into slices, like sections of an orange. 
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Defining a spherical harmonic gravitational model with a specific 
resolution at the equator provides an overly fine longitudinal resolution at 
the poles. Obviously, finer surface resolutions require higher degree models. 
The number of harmonic coefficients, however, grows as the square of the 
degree, so increasing the degree of a model is not an insignificant task. A 
more efficient model, in terms of coefficients required to define a desired 
surface resolution, would use equal size mass density patches over the surface 
of the body. Table 2.2.2-1 lists the degree and order of spherical harmonic 
models and the number of surface mass density patches required to get 
various lunar surface angular resolutions for the gravitational potential. The 
table ignores the difficult patch layout problem and resultant inefficiencies, so 
that the number of n° X n° surface mass density patches is the patch area in 
steradians (n 2 7c 2 /180 2 ) divided into the number of steradians in a sphere (4 tc). 

Table 2.2.2-1: Spherical Harmonic and Surface Mass Density Surface Resolution Comparison 


Resolution 

Spherical Harmonics 

Surface Layer 

Degrees 

Km 

Traverse Time 
(100 km Alt) 

Degree & 
Order 

# of 

Coefficients 

# of Patches 

22.5 

683 

442 s 

16X16 

285 

82 

11.25 

341 

221 S 

30X30 

957 

326 

3.297 

100 

65 s 

109 X109 

12,096 

3,795 

1.648 

50 

32 s 

218X218 

47,957 

15,190 

1.000 

33 

20 s 

360X360 

130,317 

41,253 


Table 2.2.2-1 lists the travel time to cross a patch of given size, since this 
factor is important in determining the resolution obtainable by observing a 
low-altitude orbiter. For estimation purposes, there should be at least two 
observations (or measurements) within the time it takes to fly over a given 
patch. Therefore, with 60 second Doppler count intervals of a satellite in a 100 
km altitude orbit, a spherical harmonic expansion of degree and order 30 is 
theoretically possible; more than three observations are obtained per surface 
patch. Unfortunately, this rule of thumb does not account for smaller 
spherical harmonic patches at the poles or the lack of observations during 
lunar occultations. 

As the table shows, a surface mass density model provides an economy 
in the number of coefficients estimated, and might be a preferred approach for 
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modeling or estimating the lunar gravitational field. The simulations of this 
thesis, however, were done with the spherical harmonic model since a 
surface mass density model has not been validated for PEP. Nonetheless, 
some of these other modeling techniques are discussed here for completeness. 


2.3 Localized Surface-Layer Gravitational Field Models 

Since lunar mascons act as localized gravitational disturbances, their 
modeling requires a fine degree of resolution to capture their high frequency 
content. A very high degree spherical harmonic model is required to model 
this behavior, which requires the estimation of a very large number of 
coefficients. As shown in the previous section, a more efficient method of 
modeling this behavior may be obtained through models focusing on the 
local, rather than global, behavior. In an attempt to recreate their high 
harmonic frequency behavior, mascons have been modeled by point masses, 
lens shaped mass concentrations, and gravity dipoles. In some instances, 
these mascon models have been used to model the entire lunar gravitational 
field, and in other cases, they have been combined with low degree spherical 
harmonic expansions. 


2.3.1 Point Mass Model 

The point mass model can be used to represent the behavior of an 
individual mascon. The gravitational force due to a point mass. Equation 
(2.1-1), results from taking the gradient of the following potential. 

U (f 1 ) = TZ —(2.3.1-1) 

l r i- r 2| 

Multiple mascons can be modeled by summing the potential contribution of 
each individual point mass. The resulting potential model for n-1 point 
masses becomes 


40 









Chapter Two: Gravity Field Models 


/ \ XT* Gift{ 

i=2 l r l r t 


(2.3.1-2) 


Typically, these point masses are distributed about known mascon 
locations. To recreate mascon behavior, the location, mass and depth of the 
point masses are adjusted to fit observations. Point masses are positioned 
below the lunar surface to avoid singular conditions which would result as 
the separation distance approaches zero. To globally model the moon's 
gravitational field with point masses, the model should constrain the total 
mass and center of mass of the system to known values. 


2.3.2 Surface Disk Model 


Because point mass models did not satisfactorily fit lunar orbit 
observations, scientists turned to more sophisticated surface layer 
representations. Scientists at the Aerospace Corporation [47] and the Jet 
Propulsion Laboratory [4] replaced the point mass model with a surface disk 
or lens shaped model. This model is derived from the potential of an 
ellipsoid of uniform density which is given by Equation (2.1-4) with the 
boundary of the ellipsoid used for the limits of integration. The density 
function of (2.1-4) is assumed to be constant and integrates to the mass of the 
body. The boundary condition of an ellipsoid is 



(2.3.2-1) 


where a, b, and c are the dimensions of the ellipsoid along the principal axes 
x, y and z. The surface disk or lens model uses the specialized case of an 
oblate spheroid. Figure 2.3.2-1 shows a prolate spheroid (a=b<c) and an oblate 
spheroid (a=b>c). The gravitational attraction for a point outside an oblate 
spheroid of uniform density is [34] 
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where e is the spheroid's eccentricity, so that c 2 =a 2 (l-e 2 ), and the variable k is a 
positive solution to the quadratic equation 

a 2 k 2 +(a 2 -(x 2 +y 2 + z 2 ))fc-z 2 =0 (2.3.2-3) 



Figure 23.2-1: Prolate and Oblate Spheroids 


As the thickness, c, of the spheroid approaches zero, the gravitational 
attraction of the disk is obtained. The resulting forces are then [47] 
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(2.3.2-4) 


The mass of the disk is m, a is its radius, and the x, y, and z coordinates are as 
in Figure 2.3.2-1. In the limit as the disk's radius, a, approaches zero, the disk 
shaped model approaches the point mass model. 


The Aerospace Corporation scientists used 600 surface disks of 50 km 
radii covering the lunar surface to model its gravitational field [47]. At JPL, 
they used 117 lens-shaped mass concentrations placed about 50 km below the 
lunar surface to model the gravitational field [5]. In the JPL model, the lens¬ 
shaped mass concentrations augmented the moon's central body attraction. 
Both models employed positive and negative mass disks in order to recreate 
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the observations from lunar satellites. Naturally, negative (repulsive) mass 
surface disks are not physically realizable, but they do model situations where 
a lunar surface feature is significantly less dense than its surroundings. 

Both of these surface disk models provided better correlation to lunar 
observations than previous models. Unfortunately, they did not constrain 
the lunar center of mass. When the models were converted to spherical 
harmonic models, Ji, Cn and Sn coefficients were required. Additionally, 
both models assumed a priori knowledge of the mass concentration locations 
and placed them about the moon's surface (or just below it) on a grid pattern. 
If such a priori knowledge were not used, the estimation process would 
require five terms to describe each surface disk: radius from the moon's 
center, latitude, longitude, strength (mass), and disk radius. A better model 
would have allowed the concentration's location to vary and would have 
constrained the center of mass. This, however, would have involved too 
many variables for the model to converge with the given lunar orbiter data. 


2.3.3 Gravity Dipole Model 

A gravity dipole model has also been proposed to account for the effect 
of anomalous mass concentrations upon low orbiting bodies [18]. A 
gravitational dipole consists of a mass +m separated from a fictitious mass -m 
by a distance d (Figure 2.3.3-1). In the limit as d approaches zero, m is 
assumed to get larger so that the product md remains constant. The strength 
of the dipole, D, is the result of the following limit 

n lim m(f') d = D(f ') (2.3.3-1) 

0 

Gravitational dipoles have never been shown to physically exist, but a 
distribution of these dipoles can be useful for modeling the lunar gravity 
field. 


The gravitational potential due to a gravity dipole may be written as 


U(F) = -G§ 


D(r') cos 0 
\r -f' {2 


ds 


(2.3.3-2) 
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where ds is a patch of the surface and the other elements are pictured in 
Figure 2.3.3-1. 



Gravity dipoles can be useful modeling tools because of their unique 
properties. The most useful property is that they produce a discontinuity in 
the tangential component of the gravitational field when traversing the 
dipole layer. The gravitational field normal to the dipole layer, however, 
remains continuous as the layer is crossed. This allows gravity dipoles to 
model unexplained out of plane accelerations. Using a ring of mass and 
gravity dipoles in a continuous line distribution around the lunar equator, 
Croopnick [18] successfully modeled actual disturbing accelerations beneath a 
low altitude orbiter. The dipoles were oriented normal to the equatorial 
plane, so that they were used to account for out-of-plane disturbing 
accelerations. The mass ring, meanwhile, accounted for the radial and 
tangential (in plane) components of the disturbing accelerations. 
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Orbit Propagation 


3.1 Numerical Integration Techniques 

Two classes of perturbation methods are used in celestial mechanics to 
determine precise spacecraft orbits. General perturbations generalize the 
expressions for two-body motion to include disturbing effects of other bodies 
using infinite trigonometric series expansions and integrate these series term 
by term. Special peturbations use numerical methods for deriving the 
disturbed orbit by direct integration of the rectangular coordinates or a set of 
osculating orbital elements [10]. 

Orbit propagation in PEP is performed using the numerical methods of 
special perturbations. Numerical orbit determination techniques are prefered 
because of the ease of implementation and the accuracy of solutions. The 
growing capabilities and increasing speeds of modem digital computers have 
significantly increased this method's accuracy and utility. Using Cowell's 
method, PEP integrates the equations of motion in rectangular coordinates 
fixed in inertial space. Section 3.2 discusses the units and coordinate frames 
used in this study's analyses. The equations of motion and the equations for 
the partial derivatives of motion with respect to gravity harmonic coefficients 
and initial osculating orbital elements, as used in PEP, are covered in Sections 
3.3 and 3.4. 

Two of PEP's numerical integration techniques were used for this 
study: the Adams-Moulton and Nordsieck methods. The Adams-Moulton 
constant step size integration technique with 11 th differences was used for the 
propagation of near-circular lunar orbits [15]. This technique uses predictor- 
corrector techniques to accurately extrapolate forward in time a satellite's 
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position and velocity. Its integration step size should be smaller then 1/100 th 
of the satellite's orbital period to ensure numerical stability. An even smaller 
step size must be used to accurately sample higher degree gravitational 
harmonics. A very conservative rule of thumb is to use a step size of l/(100n) 
of an orbital period to simulate an n X n degree and order spherical harmonic 
model. In PEP the step size is 2 m days, where m is a negative number. A two 
hour lunar orbit using an m = -10 step size will have over 170 steps per orbit, 
ensuring numerical stability. An m = -14 step size for the same orbit will 
have over 1,365 steps per orbit, which should adequately model a 13 X 13 
spherical harmonic gravity model. 

Elliptical orbits were propagated by the Nordsieck fifth degree variable 
step size technique [38]. Since the technique is self starting, it is used to start 
the Adams-Moulton method. This technique predicts the orbit ahead using a 
fifth degree polynomial whose coefficients are approximations to derivatives 
of the function being integrated. The integration output file therefore 
contains the satellite position, velocity, acceleration, and jerk [8]. The variable 
step size uses a smaller step near periapse, where quantities change more 
rapidly. For highly elliptical orbits, this is a more efficient integration 
technique, since constant step size methods propagate the entire orbit with 
the smallest required step size. 

Despite the integration step size used, integration quantities are written 
to an output file using a different step size. This output step size is generally 
two times the integration step size. For low frequency orbital disturbances, 
less frequent output step sizes can be used. A satellite's position and/or 
velocity are determined at specific observation times by interpolating from 
the satellite's integration output file. Everett eighth difference interpolation 
is performed on constant step size integration files. This method fits a ninth 
degree polynomial through the ten output times surrounding the 
observation time. Hermite interpolation is performed on variable step size 
integration files. This method uses a fifth degree polynomial agreeing with 
position, velocity, and acceleration at the two output times surrounding the 
observation time [8]. 
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3.2 Units and Coordinate Systems 

The numerical integration methods of special perturbations require 
that units of time and the gravitational constant are precisely defined. The 
time unit in PEP is the coordinate time day defined in terms of atomic time in 
Appendix B. The distance unit is the Astronomical Unit (AU), defined by 
setting the square root of the gravitational constant times the mass of the sun 
to the Gaussian value (see Appendix B). 

Numerical integration of the satellite equations of motion are 
performed in an inertial Cartesian coordinate system. PEP's current inertial 
coordinate system is based upon the mean equinox and equator of the earth of 
1950.0, or Julian Date 2,433,282.423. This coordinate system uses the earth's 
rotation axis at this time as the z axis, with the x axis along the 1950.0 mean 
equinox pointing towards the constellation Aries, and the y axis completing a 
right-handed coordinate system. The transformations between this 
coordinate system and those fixed in the moon and earth are discussed in 
Appendix C. 

For lunar satellite propagations, the origin of this coordinate system is 
placed at the moon's center of mass. The coordinates of perturbing bodies 
during an integration are determined from the coordinates of the earth-moon 
barycenter relative to the sun, of the moon relative to the earth, and of 
planets relative to the sun calculated by interpolation from an n-body file 
supplied from the Smithsonian Astrophysical Observatory (SAO). This n- 
body file is based upon the SAO's fit to observational data. The moon's 
coordinates are based upon formulas from Brown's lunar theory. 

To analyze the PEP propagated orbits, auxiliary software was written to 
transform the inertial integration Cartesian coordinates into a selenographic 
coordinate system. The selenographic coordinate system is also centered at 
the lunar center of mass, but its z axis is the lunar rotation axis, the x axis 
points toward the earth, and the y axis completes the right hand coordinate 
system (see Appendix C.l). One auxiliary software program computes 
selenographic orbital elements and satellite ground tracks from an inertial 
integration file. Initial conditions for lunar orbits were chosen in the 
selenographic coordinate system. The orbital angles (i, £2, co) were then 
transformed to the inertial integration coordinate frame in a second auxiliary 
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software program using the transpose of the transformation matrix used in 
the first program. 

For this thesis, inertial Cartesian coordinate frame positions will be 
referred to as X, and in selenographic Cartesian coordinates as U, where the 
two are related by the transformation matrix R(t) as follows (Appendix C.l) 

x u 

X = y U= v (3.2-1) 

z w 

Rl 1 ^12 ^13 ^1 

U = R(f)X = R 21 R 22 R 23 X= Rl X (3.2-2) 

R 3 i R 32 ^33 _ _^3 

Additionally, the selenographic Cartesian coordinates are related to the 
selenographic spherical coordinates (r,<t>,0) by the relationship 

u rcos9cos<f) 

v = rsin0cos0 (3.2-3) 

w_ rsin0 

so that the inertial coordinate frame can be related to the selenographic 
spherical coordinates used in the spherical harmonic expansion by the 
relations 

r = |x| = (X T X)^=|Q| (3.2-4) 

sin 0 = ^(xjx) cos 0 = sj\ - sin 2 0 (3.2-5) 

cose = 7zkj(R{x) sil^T^XjX) (3.2-6) 

These transformations are used in PEP's internal transformation 
routines, in the two auxiliary software programs, and in the partial 
derivatives of satellite motion equations covered in Section 3.4. 


48 





Chapter Three: Qj&il _Propapatton 


3.3 Equations of Satellite Motion 

Because of the nature of the equations of motion for a satellite, PEP 
integrates an augmented state vector. 


x; = 


A 


(3.3-1) 


where X b , is the position of body b (satellite) with respect to body 1 (the moon) 
defined by (3.2-1) and X bt is its time derivative: 

dX, 

dt 


v _ wy H/ _ r 

A w - Jj. - L 


x y z] T 


(3.3-2) 


For the following equations, X, X, and X* will all be used, although 
the reader should realize that "state" refers to the augmented state vector, X*. 

The equations of motion for a lunar satellite relative to the moon may 
be written as: 


= -GM, + H, + y ¥ + H e + other forces 

dt 2 rl 


(3.3-3) 


with the initial condition at t=to 


Y* — Y* 

A w “ A wo 


(3.3-4) 


In this thesis, the "other forces" (radiation pressure, gas leaking, 
thruster firing) are ignored, but can be included when fitting to real 
observations. In Equation (3.3-3), the Hj and H e terms are the effects of 
gravitational harmonics (zonal and tesseral) for the moon and earth 
respectively, and the ¥ term is the point mass perturbing accelerations of 
other bodies upon the spacecraft and the subscripts e, s, and p refer to the 
earth, sun and planets. 


'?'= X GM ‘ 


i=e,s,p 


X ib X, 


ib 


r 3 

r il 


(3.3-5) 


For the simulations run, the earth and sun perturbing attractions were 
included. The effect of all other bodies was considered negligible for these 
simulations. 
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3.3.1 Effect of Lunar Gravitational Harmonics 


The acceleration on a satellite due to the gravitational potential U of 
the moon is obtained by taking the gradient of the potential (2.1-3). Since 
Equation (3.3-3) already accounts for the central body attraction, the effect of 
the remaining harmonics is derived from (2.2-9) without the central body 
term. The resulting acceleration is 



where % is the ratio of p i/ibi, N z is the degree of zonals used in the spherical 
harmonic expansion and Nf is the degree and order of tesseral harmonics in 
the expansion. The recursive formulas for P n , P' n , P nm , and P' m in terms of 
the argument sin<|> which have been coded in PEP are given in Appendix A. 

The partial derivatives of the selenographic spherical angles in (3.3.1-1) 
are obtained by differentiating Equations (3.2-5) and (3.2-6). This 
differentiation results in the following relationships 


<?sin0 

3X„ 


= JL - ^-sin 0 

r bl 


(3.3.1-2) 


-sin 6 


R x cos 6 f cos 0 

r cos 0 cos 0 1 dX bl 


-2-cosd (3.3.1-3) 

rli 


cos 9 


de y 

dX blJ 


T 

R 2 sinfl r ^ cos 0 ^ _^i± sin Q 

rcos0 cos 0 dX bl J r\ x 


(3.3.1-4) 
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Equations (3.3.1-3) and (3.3.1-4) can be simplified by multiplying the 
first by -sin0, the second by cos0, and then adding the two equations. This 
yields the following formula for the partial derivative with respect to 0. 


de ] 

T 

_ 1 | 

Ux w J 

r bl cos 0 1 


(3.3.1-5) 


Equation (3.3.1-1) is singular at $ = ±90° for PEP's algorithms. This is 
generally not a problem since it is highly unlikely that a satellite numerical 
integration step would land exactly over a pole. One simulation, however, 
used a polar location as an initial condition and PEP was unable to propagate 
its orbit. The initial mean anomaly for this orbit was altered by one degree 
and the integration proceeded without any further difficulties. 


3.3.2 Effect of Other Gravitational Body Harmonics 

PEP has the capability to include other (non-central) gravitational body 
perturbing effects upon a satellite's motion. This enables PEP to handle 
spacecraft fly-by missions. PEP can include the effects of a destination, or 
target, body's harmonics upon the motion of a satellite traveling towards one 
body but which is within the sphere of influence of another body. Using this 
feature, PEP can also include the earth's J 2 harmonic upon a lunar orbiting 
satellite. Higher order terms can also be included for integration accuracy, but 
are rarely required. 

The effect of earth perturbations on the motion of a lunar satellite 
relative to the moon (H e from (3.3-3)) is the difference between the effect of 
the earth on the satellite (H*, e ) and the effect of the earth on the moon (H^). 
H be is given by Equation (3.3.1-1) with a change of subscripts, replacing 
subscript 1 with e. The effect of earth harmonics upon the moon is calculated 
by replacing the GMj terms in (3.3.1-1) by GMc, (the subscript c refers to the 
earth-moon barycenter), the subscript 1 by e, and the subscript b by 1. The H e 
term in (3.3-3) is then given by 

R, - Rbe - H„ (3.3.2-1) 
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Since low altitude lunar orbits and eccentric orbits with apolune on the 
far side of the moon were simulated, no earth perturbing effects were 
included. When processing real observations to estimate the lunar 
gravitational field, earth's J 2 harmonic should be considered since it is three 
orders of magnitude larger than the other earth harmonics. 


3.3.3 Effect of Mass Concentrations 

For this thesis, PEP was modified to include the effects of mascons 
upon lunar orbiting satellites. Since these modifications have not been 
integrated into the SAO's version of the software, the Draper Laboratory 
modified version will be referred to as PEP-D. PEP-D implements mascons by 
modeling them as point masses (Section 2.3.1). 

Due to mascon model implementation, the equations of motion for a 
lunar satellite relative to the moon (3.3-3) become, in PEP-D, 

= -GM l (1 -v)^ L + H l +^ , + H e +K (3.3.3-1) 
dt T bl 

with the same initial conditions, (3.3-4). In Equation (3.3.3-1), the K term is 
the effect of all of the mascons and u is their total mass fraction. 

The mass of each mascon and its selenographic position in spherical 
coordinates is entered into PEP-D as the program is initialized. The mass or 
strength of a mascon is input as a fraction of the total central body mass and 
both positive and negative values are allowed. If N* is the total number of 
mascons, the total mass fraction of the mascons is then calculated by 

N k 

v = '£ i m i 0< U<1 (3.3.3-2) 

i=l 

The original central body term from (3.3-3) is reduced by the factor (1-u) 
to conserve mass in the lunar system. This will allow the lunar mascon 
model to behave identically to the central force model far from the moon. 
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The acceleration on a satellite due to the gravitational potential of N* 
mascons is obtained by taking the gradient of the point mass potential (2.3.1- 
2). The resulting acceleration is 


N k 


K = -GM { 


7W;X, 


bi 


i=l 


bi 


(3.33-3) 


From the input spherical mascon selenographic coordinates, PEP-D 
determines their inertial Cartesian coordinates X )7 using Equation (3.2-3) and 
the inverse of Equation (3.2-2). The coordinates of the lunar satellite relative 
to mascon i and their separation are then 



= X W -X fl 


(3.3.3-4) 


*bi 



(3.3.3-5) 


Unfortunately the mascon implementation in PEP-D does not 
constrain the lunar center of mass. As mentioned in Section 2.3.2, when 
converting a surface layer model with an unconstrained center of mass to a 
spherical harmonic expansion, first degree harmonic coefficients ]\, Cn, and 
Sn need to be determined. Rather than modify PEP-D to estimate first degree 
coefficients, this thesis uses mascon models in which the lunar center of mass 
is not disturbed. If any mascons are used, PEP-D calculates the central body 
center of mass by the formula 

N k 

x e .„. = (3-33-6) 

1=1 

When the somewhat arbitrary mascon "truth" model was created (Section 
6.4), this feature was used to determine the placement and mass of 
"balancing" mascons to preserve the lunar center of mass. 


3.4 Partial Derivatives of the Satellite Motion Differential Equation 

The partial derivative of Equation (3.3-3) with respect to a parameter a 
yields the variational equations that are numerically integrated in PEP along 
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with the equations of motion. The parameter a could be one of the orbit's 
initial conditions, the harmonic coefficients of the gravitational field, or other 
parameters of interest. Taking into account the effects of mascons, the partial 
derivative of Equation (3.3.3-1) with respect to the same parameter a yields 
PEP-D's variational equations 


( A 
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da 
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bl 
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dK 
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(3.4-1) 


with the initial condition at t = to 


dX bl _ dXw) 
da da 


(3.4-2) 


Unless estimating the mass of the moon, the first term in (3.4-1) is zero. 
Additionally, since the mascon masses were not estimated in this thesis, the 
ratio i) was considered constant. 


The effect of perturbing body attractions (point mass approximation) 
upon the partial derivative with respect to a is obtained by differentiating 
Equation (3.3-5), yielding 
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(3.4-3) 


3.4.1 Effect of Lunar Spherical Harmonics on the Partial Derivatives 

The effect of the lunar gravitational harmonics upon the partial 
derivatives of satellite motion with respect to a parameter a is obtained by 
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differentiating Equation (3.3.1-1) term by term. After differentiating, like 
series summations can be regrouped for numerical algorithms. The resulting 
equation takes the form 


9H, a, 3 ( GM t) , CM, <= f 1 
da GM, da r 2 u 1 J 


(3.4.1-1) 


where Z is the series expansion for the zonal terms and f is the series 
expansion for the tesseral terms. The first term in (3.4.1-1) is set to zero unless 
attempting to estimate the mass of the moon. Z, the zonal series expansion, 
is provided by the following expression: 
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(3.4.1-2) 


The first term in the above equation is zero since the zonal coefficients are 
constant. 


After grouping like terms, the series expansion for the tesserals, T , 
reduces to the following expression: 
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(3.4.1-3) 


The first two terms in this expansion are zero since the tesseral spherical 
harmonic coefficients are constant. The expansion for Ti and T 2 are then 
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(3.4.1-4) 
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nm 


The recursive formulas for calculating P n' P n r Pn> Pnm ' ^nm' an< ^ T 
in terms of their argument sin<]) are included in Appendix A. 

Additionally, based upon the transformation between selenographic 
spherical coordinates and the inertial frame (3.2-4), (3.2-5) and (3.2-6), the 
following partial derivatives are obtained. 
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PEP has the option of using different degrees of spherical harmonic 
models for calculating the partial derivatives in Equation (3.4-1) than for 
calculating the satellite motion in Equation (3.3-3) or (3.3.3-1). Since 
calculating the effect of the harmonics is such a tedious process, calculating 
the partial derivative's harmonic effect with a lower degree spherical 
harmonic model can save quite a bit of computer time. This option works 
well for high altitude earth satellites where J 2 dominates the other harmonic 
effects. This option was found impractical for low altitude lunar satellites, so 
the same degrees N z and N t were used as the summation upper limits in 
Equations (3.4-1) as in (3.3.3-1). 


3.4.2 Effect of Earth Spherical Harmonics on the Partial Derivatives 

The effect of earth harmonics on the partial derivatives of satellite 
motion with respect to a parameter a is calculated by differentiating Equation 
(3.3.2-1). This calculation involves differentiating two factors of the form of 
Equation (3.4.1-1) since the effect of the earth harmonics on the moon must be 
subtracted from the effect of the earth harmonics upon the lunar satellite. For 
this thesis, the effect of earth harmonics upon the partial derivatives were 
neglected in Equation (3.4-1) because all of the orbits propagated were low 
altitude ones: 



(3.4.2-1) 
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3.4.3 Effect of Mass Concentrations upon the Partial Derivatives 


The effect of mass concentrations on the partial derivatives of satellite 
motion with respect to a parameter a is obtained by differentiating Equation 
(3.33-3) assuming that the total lunar mass remains constant. 
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(3.4.3-1) 


The above equation was not coded into PEP-D. The partial equations 
were calculated on orbit fitting runs and all of the orbit fitting was done to 
spherical harmonic models. The spherical harmonic plus mascon truth 
model was run to generate the satellite observations for each gravitational 
sensing method, for which partial derivatives were not required. If PEP-D 
was used to estimate coefficients in a spherical harmonic plus mascon model, 
this equation would be required and could be used in an estimation fit which 
varied the mascon strength and location. 


3.4.4 Effects of Initial Conditions on the Partial Derivatives 


The partial derivatives of the satellite's initial conditions with respect 
to a from Equation (3.3-4) (used as initial conditions (3.4-2)) are zero unless 
the parameter itself is a satellite initial condition. For the cases where a is a 
Cartesian coordinate initial condition, the partial derivative is obtained from 
the following relation. 


Mm = j 


1 0 0 0 0 0 
0 1 0 0 0 0 
0 0 1 0 0 0 
0 0 0 1 0 0 
0 0 0 0 1 0 
0 0 0 0 0 1 


(3.4.4-1) 


Rather than solve the partial derivatives with respect to inertial 
Cartesian initial conditions, PEP uses the partial derivatives of the initial 
osculating elliptical orbital elements (ao, eo, io/^o, «*)/ Mo) or (ao, eo, io/ ^o. 
(Q+0))o, (Q+co+M) 0 ). The osculating orbit is the one which would result if the 
disturbing forces in Equations (3.3-3) or (3.3.3-1) were instantly "turned off 
and the satellite's motion continued along the two-body orbit defined by the 
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central body force. These orbital elements, like the Cartesian coordinate state 
vector, completely describe the satellite's initial state. Calculating the partial 
derivatives with respect to the initial osculating elliptical orbital elements is 
numerically more efficient because only the partial derivative with respect to 
the initial osculating semi-major axis grows with time. Reference [6] contains 
the formulas for the partial derivatives of the initial conditions with respect 
to the initial osculating elliptic orbital elements in the integration coordinate 
frame. 


3.4.5 Checking Partial Derivatives by the Difference Method 


To verify that partial derivatives were being calculated correctly, they 
were checked using a finite difference method. This method was used to 
check both position, X, and observation partials. For these checks, orbits 
and/or observations with partials were generated with two different values of 
the parameter a (ft] and 012 ). The finite difference method then verifies the 
partial derivatives of the quantity of interest, T (scalar position coordinate or 
theoretical observation), by verifying the following equality. 


dr 


dr 


da 

«=i 

+ da 

a=«2 




at-a 2 


(3.4.5-1) 


This check allows coding errors to be detected and verifies the proper 
execution of the program. Once coding errors are fixed and the program is 
being run properly (appropriate control inputs are used), the orbit fitting 
process generally converges to a solution. If the program is running properly 
but has difficulty converging upon a solution, this indicates that parameters 
are too highly correlated, the observations do not provide adequate 
observability, or the fit model is an inadequate representation of the observed 
behavior. 
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Satellite Observations 


4.1 Observation Simulation and Processing 

Knowledge of the world in which we live is based upon observations, 
or measurements, of physical phenomenon. Theories that attempt to explain 
and model natural occurrences are developed to explain previous 
observations and are validated by accurately predicting future behavior. The 
discovery of the planet Neptune is an example of the impact observations 
have upon the development and refinement of theoretical models. After the 
discovery of Uranus, astronomers could not reconcile observations of the 
planet's motion with theoretical predictions. English astronomer John Couch 
Adams and French astronomer Urbain-Jean-Joseph Le Verrier, independently 
studying the motion, concluded that Uranus' behavior was due to a planet 
beyond it. Using a new model, both astronomers were able to predict the 
location of this ultra-Uranian planet, and shortly thereafter Neptune was 
discovered [10, pp. 472-473]. 

The Planetary Ephemeris Program (PEP) was designed to process 
observations of the sun, moon, planets, stars, and spacecraft and to further 
scientific knowledge through this observation processing [8]. PEP processes 
observations of interplanetary, earth, and lunar spacecraft using observation 
files. These files, based on actual astronomical observations or simulated 
observations, are made up of observation series. Observation series are 
limited to a single sensing method (Section 4.2) and a pair of observation 
types (Section 4.3). Multiple observation methods, types, periods, and 
frequencies are accommodated within PEP by either including more than one 
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observation series in an observation file or by using more than one 
observation file at a time. 

To optimally estimate a spherical harmonic expansion fit to the lunar 
gravity field for this thesis, PEP generated simulated "truth" model 
observations based upon the observation geometry and measurement type 
being evaluated. Observation series for each observation geometry and 
measurement type were created using PEP's dummy observation feature 
where output observation files were created from knowledge of earth-based 
observing sites, and observing and observed satellite states determined from 
numerical orbit propagations which used the gravitational "truth" model. 
These observation series included each method's associated measurement 
error - an indication of the statistical accuracy of any given measurement. 
The error values used were based upon published or calculated instrument 
capabilities. During the fitting process, the observation files were read and 
these errors were used to weight the observations. 

Treating the "truth" model observations as real, PEP fit a spherical 
harmonic model to these observations. Using first guesses for the fit model 
coefficients and satellite initial osculating orbital elements, the equations of 
motion and the equations for the partial derivatives of the motion with 
respect to these quantities were numerically integrated. From the numerical 
integrations, the theoretical observation values were determined. 
Observation residuals, the difference between the "true" observations and the 
theoretical values calculated, and the observation partial derivatives were 
computed and written to an observation output file. This file was later used 
to update the guess for the gravitational coefficients and orbital initial 
conditions. 

There was no measurement noise added to the "truth" model 
observations in this thesis. White measurement noise could be added using a 
random number generator, with any non-zero measurement biases estimated 
by PEP. 

Both "truth" model and theoretical observation values depend upon 
the location of the observing site, the observation time (signal receive time), 
the location of the observed body at reflection time, and any distortion factors. 
PEP determines the receiving site's position or state at signal receive time in 
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the inertial Cartesian coordinate frame based upon integration or n-body files. 
For observation processing in cis-lunar space, position is determined relative 
to the center of the earth. For extremely precise observations, PEP can 
perform the light time iterations with the center of mass of the solar system 
as the center of the inertial coordinate system. 

The distance between the observed body at a guessed reflection, 
transpond, or transmission time and the observing site at receive time is 
converted to light time, and the reflection, transpond, or transmission time is 
adjusted until the difference between the reflection, transpond, or 
transmission time and receive time is sufficiently close to the distance 
between the bodies converted to light seconds. 

For two way observation signals, this process is performed twice. The 
first iteration is used to determine the reflection time and state of the 
observed body. The second step uses the reflection time (less any systemic 
delays) and state of the observed body to compute the signal source's send 
time and position. Bent pipe observations perform these light time iterations 
as many times as necessary to recreate the path of the electromagnetic signals 
from signal receive time back to its origin. 


4.2 Observation Methods 

Astronomical observations are collected using several different 
methods. The most common observations are taken from earth sites. These 
can be observations of celestial bodies or man-made spacecraft in earth, lunar, 
or interplanetary orbits. Additionally, satellite based observations are made of 
other satellites or of sites on the earth or lunar surface. Finally, bent pipe 
observations are made in which a signal is passed among multiple satellites 
and sites. PEP can simulate observations using each of these methods and can 
also model physical effects which affect these observations. 
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4.2.1 Earth-Based Observations 

PEP contains the coordinates of several earth based sites to process 
astronomical observations. DSN sites are included in cylindrical coordinates 
and most other locations are in spherical coordinates. Observation sites not 
already included can be added to PEP. This will permit the creation of lunar 
observation sites. If a lunar base or navigation site were established, PEP 
could process lunar-based observations in the same manner used to process 
earth-based observations, after some straightforward software modifications. 

When an earth-based site is involved in an observation, the site's 
inertial location with respect to the center of the earth must be determined to 
calculate the theoretical value of the observations. Given an observation 
receive or send time in UTC time (as disseminated by the U.S. Naval 
Observatory Time Service Radio Station, WWV, and other time services), 
PEP determines the coordinate time (CT), UT1 time, and earth wobble 
coordinates from look-up tables. International Astronomical Union (IAU) 
expressions for the sidereal time as a function of UT1 time and the earth 
precession-nutation matrices as functions of CT time are then evaluated. PEP 
transforms the observing site's earth fixed coordinates to the integration 
frame coordinate system using the transformations described in Appendix 
C.2. PEP can also calculate the partial derivatives of the observing site's 
integration frame coordinates with respect to cylindrical or spherical 
coordinates. This feature allows PEP to improve the estimate of the 
observation site's location in earth fixed coordinates as part of the process. As 
lunar sites are established, PEP can therefore survey the site's location by 
processing satellite-, earth-, and lunar-based observations. 


4.2.2 Satellite-Based Observations 

PEP also allows observations when a satellite is the signal receiving 
and/or sending site. In this case, PEP uses the satellite's integration file to 
determine its state and the partial derivatives of the state at receive, 
bounce/transpond, or send time. When using a lunar satellite for 
observations, the satellite's state and partial derivatives are translated from 
lunar- to earth-centered coordinates. The lunar satellite's state and any partial 
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derivatives at receive, bounce/transpond, or send time are determined using 
the interpolation methods described in Section 3.1. 


4.2.3 Bent Pipe Observations 

PEP also processes observations in which signals pass among several 
satellites and sites. This coding has been used to recreate two different 
ranging observations [8]. In the first method, an earth site transmitted a 
ranging signal to an earth satellite which then transponded the signal to a 
second earth satellite. This second satellite then transponded the signal to an 
earth receiving site. The second bent pipe method used the same signal path 
from the originating earth site through the second satellite, but then the 
signal was transponded back to the first satellite and then transponded to an 
earth receiving site. In both cases the earth sending and receiving sites can be 
the same or different sites. PEP-D can process these same bent pipe 
observations methods for lunar orbiting spacecraft with earth-based sending 
and receiving sites. Unfortunately these bent pipe methods currently only 
process ranging observations (Section 4.3.6). To accommodate bent pipe range 
rate measurements, PEP will have to be modified. 

A bent pipe method using two two-way coherent links has been 
proposed to observe satellite motion on the far side of the moon (Section 1.4). 
For this method, an earth-based observation site sends a high frequency signal 
to a lunar satellite. This Doppler shifted received frequency is used to 
generate a medium frequency signal for a two-way coherent Doppler loop 
between two lunar satellites. The returned Doppler shifted medium 
frequency is then modulated onto the received high frequency signal and 
transmitted to earth, completing a second coherent loop [40]. Since this 
observation method was not available for range rate measurements, it was 
simulated by using two separate, unrelated coherent Doppler loops. Since the 
proposed method depends on a single frequency source, the simulated 
method used a perfect, non-drifting frequency source for the satellite-to- 
satellite loop. This method obtained range rate observables equivalent to the 
proposed method. 
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4.2.4 Observation Interruptions 

All observations involve the transmission of an electromagnetic signal 
(light or radio waves) between sending and receiving sites and satellites. 
Observations therefore depend upon a clear line-of-sight between the sites 
and satellites. Multiple signal path methods depend upon a clear line-of-sight 
for each path. For earth observations of a lunar satellite, the line-of-sight can 
be interrupted by either the earth or the moon. The first case occurs if the 
line-of-sight passes below the horizon. In the second case, the moon occults 
the line-of-sight when the lunar satellite passes behind the moon. PEP 
models both of these observation interruptions and deletes observations from 
dummy "truth" model observation series. 

Given the radius of the occulting body, PEP determines whether the 
satellite is occulted by its central body using vector descriptions of the 
observing site, observed satellite, and central body locations. When the 
observed body becomes occulted by its central body, observations cease until 
the observed body returns to view. For the gravitational sensing methods 
simulated, this feature was used with a lunar radius slightly larger than the 
lunar radius to eliminate poor quality observations and account for satellite 
acquisition difficulties at the edges. 

Given a limiting elevation angle, as depicted in Figure 4.3.1-2, PEP 
determines whether the line-of-sight has passed below the observing site's 
horizon using vector descriptions of the line-of-sight and the vector normal 
to the observing site. PEP discontinues observations while the line-of-sight is 
below the required elevation angle. This feature was not used in the lunar 
gravity model estimation simulations since there is at least one observing site 
facing the moon at all times (these sites are spaced around the earth's 
longitude). For this thesis, the simulations used a single observing site which 
made observations through a "transparent" earth rather than using several 
observing sites and simulating the handing off of observation responsibilities 
from site to site. 
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4.3 Observation Types 

In the past, astronomical measurements were limited to optical 
sightings of celestial bodies involving angular measurements and sighting 
times. Angular sighting observation pairs include azimuth and elevation, 
right ascension and declination, and meridian crossing and elevation angle. 
More recent astronomical measurements involve electromagnetic waves 
being transmitted through space. Electromagnetic signals may be reflected, 
transponded, or transmitted by man-made bodies. In addition, natural bodies 
are also used to reflect radar signals. These more recent observation 
techniques can provide precise interferometric angular information as well as 
accurate range and range rate measurements. Independently of the 
observation method used, PEP determines the signal path from send to 
receive time using the light time iterations described in Section 4.1. PEP 
models any effects which could bend or distort this signal path based upon the 
type of observation. The aberration of light, the Doppler shift in frequency, 
atmospheric refraction, ionospheric distortions, general relativity, and 
interplanetary plasma effects are all modeled in PEP to determine the proper 
signal path. The adjusted signal path between sending and receiving sites, q, 
is then used to recreate observations. 


4.3.1 Azimuth-Elevation Observations 

Early astronomical sightings measured the angles between the apparent 
line-of-sight to the target body and a reference frame. Azimuth and elevation 
angle observations at the observing site use the vector normal to the 
observing site and the plane tangent to it to describe the location of the 
observed body. The elevation angle is the angle between the line-of-sight and 
its projection in the tangent plane. The azimuth angle is defined as the angle 
in the tangent plane between north and the line-of-sight's projection. These 
vectors and angles are depicted in Figures 4.3.1-1 and 4.3.1-2. 

In earth-fixed Cartesian coordinates, the vector p is [0 0 1] T and the unit 
normal, n, is defined by the longitude 0 and geodetic latitude <{) as 
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tl = 


COS 0 COS <j) 

sin 0 cos 0 
sin0 


(4.3.1-1) 


PEP transforms these vectors to the inertial coordinate frame in which q is 
calculated using the transformations of Appendix C.2. 



Figure 4.3.1-1: Azimuth-Elevation Angle Vectors in the Meridian Plane 

The vector m, which points along the meridian from the observation site 
towards the north, serves as the north reference direction. This vector is 
defined by the following equation: 

m - unit\p - (p • n)n] (4.3.1-2) 



Figure 4.3.1-2: Azimuth-Elevation Angles and the Tangent Plane 


The projection of the signal path, q, in the observation site's tangent plane, 
q p/ is used to define the azimuth angle and is obtained by the equation 
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(4.3.1-3) 


Additionally, Figure 4.3.1-1 shows the projection of the signal path in the 
meridian plane, q m . 


The following equations are used to calculate the observed azimuth 
and elevation angles. Two-argument arctangent routines should be used for 
improved accuracy and to distinguish the proper quadrants [8]. 


elevation = sin 




(4.3.1-4) 


azimuth = tan 



(4.3.1-5) 


4.3.2 Right Ascension-Declination Observations 


When an object is observed against a star background, right ascension 
and declination angles can be determined from the object's relationship to 
catalogued stars. Right ascension and declination are angles referred to an 
inertial Cartesian coordinate system centered within the observing body or 
the true equinox and equator of date. If the line-of-sight vector q is given in 
the relevant Cartesian coordinate system, [x^ y fl z^] T , the right ascension and 
declination angles are calculated from the following relations, where once 
again two-argument arctangent routines should be used for improved 
accuracy and quadrant determination [8]. 


declination = sin 




\Yi\J 


(4.3.2-1) 


right ascension = tan 


-i 


f y± 

K X <U 


(4.3.2-2) 
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4.3.3 Meridian Crossing-Elevation Angle Observations 

Some of the earliest observations recorded the time and elevation of a 
celestial body's passage across the earth's meridian. These observations were 
performed by constraining the sighting instrument in the plane of the 
meridian and measuring the elevation angle as the observed body crossed the 
plane. Specifying the time of meridian crossing uses the rotating earth to 
determine one of the angular components of the line-of-sight from the 
observing body. 



Figure 4.3.3-1: Elevation Angle at Meridian Crossing 

This method is similar to the azimuth-elevation observation with the 
elevation angle constrained to the meridian plane. It is also related to the 
right ascension-declination observation because the local sidereal time at 
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meridian crossing, converted to radians, is the right ascension. Based upon 
the observing body's orbit and rotation, meridian crossing time and elevation 
angles are converted to right ascension and declination angles [8]. 


4.3.4 Satellite Look Angle Observations 


When a satellite is the observing site, PEP defines the azimuth and 
elevation angles in the satellite's pitch-roll-yaw coordinate system. The yaw 
axis, y, points from the satellite to the center of the body which it orbits. The 
roll axis, f, lies in the orbital plane normal to the yaw axis, making an acute 
angle with the satellite's velocity vector. The pitch axis, p, is normal to the 
orbital plane and completes the right hand coordinate system. This 
coordinate frame is defined for a lunar satellite's position and velocity vectors 
by the relations 


y = -unit X bl j, p = unit 


*bi x Xbi 


SK ✓V Ak . _ v 

r = yxp (4.3.4-1) 



From this coordinate system, PEP's azimuth angle is defined as the 
angle about the pitch axis and its elevation angle is the rotation angle from 
this point to the line-of-sight. These angles are depicted in Figure 4.3.4-1 and 
the formulas for PEP's azimuth and elevation angles are given below [8]. 


elevation = sin 





(4.3.4-2) 
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azimuth = tan 


q *r ' 


q*y) 


(4.3.4-3) 


4.3.5 Interferometry Observations 


PEP can also determine angular measurements from interferometer 
measurements. The angle between the incoming wave and the line 
connecting the two observing sites can be determined when two sites receive 
an electromagnetic signal and determine its arrival time difference. 



Figure 4.35-1: Long Baseline Interferometry Measurement Geometry 


If the signal source (such as a radio star) is sufficiently far from two 
receiving sites separated by a distance d, then the two signal paths are 
considered parallel. If the signal arrives at site 2 at a time At after it reaches 
site 1, then this signal has traveled the additional distance cAt, where c is the 
speed of light. The angle y is then determined from the relationship 


y/ = cos 



(4.3.5-1) 


The theoretical value of an interferometer measurement of a satellite 
is the difference between the light times from the two observing sites at the 
same receive time and the satellite at the signal transmission times. 
Detecting smaller At's and increasing the separation between the observing 
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sites to even intercontinental distances increases the accuracy of 
interferometric observations. The use of hydrogen maser atomic clocks 
allows angular measurements with milli-arc second accuracy. Previously 
discussed angular measurements are limited to approximately an arc second 
of accuracy. PEP is coded to simulate and process these observations and 
model the bias between the clocks at the two receiving sites [8]. 


4.3.6 Range Observations 

The range between two bodies can be determined from the time delay 
of a signal sent between the bodies. Range measurements depend upon a 
precise knowledge of when the signal was sent and when the signal was 
received. Since it is very difficult to synchronize two separated clocks, one¬ 
way ranging measurements are not often used, except for multi-GPS satellite 
observations where receiving site clock error is measured. A two-way signal 
provides a more accurate single satellite range measurement since a single 
clock is used to measure the time delay. 

Two-way range observations can be obtained by sending a pulsed 
electromagnetic signal between two bodies. The time it takes the pulse to 
return to the sending site, less any known delays, reveals the range between 
the bodies calculated rigorously from light time iterations. Since the 
electromagnetic signal travels at the speed of light, the range is calculated 
from the compensated time delay. At, as 

R = — (4.3.6-1) 

2 

This two-way range observation can also be obtained with the 
transmission of a continuous sine wave signal by shifting the phase of the 
sine wave 180° at a certain repetition interval according to a coded pattern. 
Since this phase shifting is preserved as the signal passes through space, 
transponder electronics, or is reflected off of an observed body, the receiving 
site can recreate the phase shifting coded pattern from the received signal. By 
correlating the sent coded pattern with the received pattern over time, the 
receiving site can determine the signals round trip time delay. Using this 
phase shift keying technique on a continuous sine wave signal allows the 
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signal to be used for both range and range rate observations, as well as for 
communications, if desired [45]. 



Two-way time delay measurements in PEP are obtained directly from 
knowing the site and body positions at the appropriate times plus any 
transponder delays. Transponder time delays can be determined prior to 
launch or can be estimated from the measurements. Typically both methods 
are used by separating the time delay into known (previously measured) and 
unknown (to be estimated) parts. 


4.3.7 Range Rate Observations 

Range rate observations measure the rate at which the observation site 
and observed body are approaching one another. These observations take 
advantage of the Doppler effect upon electromagnetic signals. Since this 
observation measures the change in frequency between send and receive 
time, it has the same difficulty with one-way observations as range 
measurements. Coherent two-way signal paths, however, provide very 
accurate measurements of the range rate between two bodies. 

The Doppler effect is a shift in the frequency of an electromagnetic 
wave radiated, reflected or received by an object in motion. This frequency 
shift is a result of the expansion or compression of electromagnetic waves 
along the direction of a moving source [43]. This compression in the 
direction of motion is visualized in Figure 4.3.7-1 below. For astronomical 
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range rate measurements, the frequency of the signal sent from site 1 to site 2 
is shifted due to the relative velocity of site 1 at send time and site 2 at receive 
time along the line-of-sight. From the light time iteration process mentioned 
in Section 4.1, PEP determines this relative velocity from the state of the two 
sites at the two times in an inertial reference frame. 



Figure 4.3.7-1: Waves Radiated from a Stationary and Moving Source 


For a range rate of R, the one-way Doppler shift of a frequency source, 
f s , is given by the formula 

A/ = j (4.3.7-1) 


R is positive when the bodies are moving apart and the signal's frequency is 
decreasing due to the Doppler shift. This formula is used to determine the 
frequency received or reflected at the observed body. For radio wave 
transponder observations, the received frequency is multiplied by a rational 
factor q, so that the body transponds a different frequency, f S 2 , from the signal 
it receives. 


f s2 = if, 


( 


1 -* 




(4.3.7-2) 


Once again, after transmission, this frequency shifts due to the Doppler effect 
by (4.3.7-1) where f s is the send frequency from (4.3.7-2) and the range rate 
depends upon the new send and receive time positions and velocities. 

For approximately instantaneous send and receive times with no 
transponder frequency shifting (q = 1), the round trip Doppler shift is 
approximated by the single expression [43] 
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Af = -2/ s — (4.37-3) 

C 

The NASA Deep Space Network (DSN) measures a coherent count of 
the zero crossings of the Doppler signal, refered to as Integrated Doppler. A 
DSN site sends a continuous sine wave signal towards a spacecraft at a 
frequency controlled by an atomic frequency standard. The receiving DSN 
station can be the same or different from the sending site because of the 
network's precise synchronization. The receiving station pre-multiplies the 
transmitted frequency by the spacecraft's transponder translation factor and 
then subtracts this frequency from the received signal. Any system biases are 
then added to this differenced frequency, if known, or they can be estimated 
in the fitting process. A counter at the receiving site is incremented for every 
positive traveling zero crossing of the differenced frequency. This counter is 
read at uniform intervals, At, to determine the interval's Doppler frequency 
shift in cycles per second (Hz). Additionally, a resolver is used to give 
fractional cycle resolution of the zero crossings [36]. 

Within PEP, the theoretical value of this observable is the difference in 
round trip phase delays at the count starting and ending receive times 
multiplied by the sending frequency. PEP also accounts for any transponder 
frequency translation when calculating the theoretical value of the 
observation. The exact formula for the Doppler count observable is coded in 
PEP [8] and is approximately equal to the instantaneous Doppler shift over the 
count interval. 


4.4 Partial Derivatives of Satellite Observations 

The partial derivatives of the theoretical value of an observation are 
required, along with the observation residuals, to estimate orbital initial 
conditions and parameters. The theoretical value of an observation is a 
function of the receiving site coordinates at signal receive time, the observed 
body's coordinates at transpond time, the sending site's coordinates at send 
time, and other parameters. If bent pipe observations are used, the theoretical 
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value is also a function of any additional body 7 s coordinates at the appropriate 
times. For a given observation, T, this functional relationship is denoted by 

r(t r )=f(F r jt r )X,M,),x;jt s ),... ,p) <4.4-i) 

where t r is the receive time, r.s. is the receiving site, t.s. is the transponding or 
reflection site, t t is the transpond or reflection time, s.s. is the sending site, and 
t s is the send time. The ellipses denote that T may be a function of other 
coordinates at other times, depending upon the observation method. P is the 
vector of parameters besides motion which affect T. 

As the theoretical observations are determined, their partial 
derivatives with respect to the coordinates of interest can also be calculated. 
Depending upon the type of observation, PEP also calculates the partial 
derivatives of the observation with respect to the parameters, P , which affect 
the observation, such as measurement biases and transponder delays. The 
partial derivative of the theoretical observation, T, with respect to a parameter 
of interest, a, is then calculated by the chain rule. 

3T d£ <?X* S JT dr <?x;, , + ^lM (44 . 2) 

da dX', da dX’ ls da dX"„ da dp da 

The parameter a would be any parameter to be estimated, such as orbital 
initial conditions, gravitational harmonic coefficients, observing site 
coordinates or observation biases. The partial derivative of a satellite's 
coordinates are determined by interpolation from the satellite's integration 
file using the methods discussed in Section 3.1. The partial derivative of a 
site coordinate is a function of the site's spherical or cylindrical coordinates. 

Once calculated, these partial derivatives are written to an observation 
output file for each receive time for each observation type in an observation 
series. The theoretical observations, their partial derivatives, and the 
observation residuals for each observation series on the output file can then 
be used to calculate adjustments to the parameters a. 
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Chapter Five 


Parameter Estimation 
and 

Prediction of Uncertainty 


5.1 Parameter Estimation 

Estimation is the process of extracting information from a collection of 
observations (measurements) to develop a better understanding of the 
observed behavior. For this thesis, the parameter estimation process attempts 
to obtain the best set of lunar gravitational harmonic coefficients. The values 
obtained in the estimation process are best in the sense that the estimated 
model's generated observations provide the closest match with the truth 
model simulated observations. This does not mean that a better fit could not 
be realized with more or different "truth" model observations. It also does 
not imply that a better fit could not be obtained with a better model of the 
behavior. The estimation process additionally provides insight into the 
uncertainty of the estimates' ability to model the behavior. This information 
can be used to analyze and evaluate the estimated model. 


5.1.1 Observation Vector as a Random Variable 

Using the gravitational "truth" model, observations were generated for 
several proposed sensing schemes. These observations are dependent on the 
lunar gravitational field and orbital states of the satellites. Together these 
observations form the N-dimensional vector z where N is the total number 
of measurements. 
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With initial guesses for 1) the lunar gravitational field coefficients, 
2) satellite initial conditions (initial osculating orbital elements for each 
satellite used in the sensing scheme), and 3) measurement biases, theoretical 
values of the estimated model observations, the difference between the two 
sets of observations, and the partial derivatives of the estimated model 
observations are generated. The estimated model observations' theoretical 
values are a function of the observation time, the estimated model 
parameters, and other factors, such as the positions of the observing sites. 
These other factors are treated as known quantities for the estimation process. 
Let dig be the vector of unknown gravitational parameters and ft be the 
vector of unknown measurement biases and any other unknown parameters. 
For a sensing scheme with N& satellites, there are N& vectors a n of unknown 
orbital initial conditions. If a is the vector of all of these unknown system 
parameters, then the theoretical observation at time t is a function of t and a. 


[h - Jk. 

|C 21 
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... 

w 

(5.1.1-1) 
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(5.1.1-3) 


Due to imperfect knowledge of the gravitational field coefficients and 
satellite initial conditions and the difference between the fit and truth 
models, the theoretical observations do not perfectly recreate the "truth" 
model observations. For each measurement time, t, there is a measurement 
error, 0(t), between the "truth" model and theoretical observations. With the 
"true" values for the parameters a, the remaining measurement error, 0, is 
assumed to be a zero-mean random variable (disregarding the difference 
between the fit and truth models). Since 0 is the result of several 
independent random causes, it will be assumed to be normally distributed by 
the central limit theorem. The collection of N measurement errors over the 
observation period is then 6. 

From the previous definitions, it follows that the vector of 
measurements is a random vector composed of a deterministic value, the 
theoretical observation, and a random quantity, the measurement error. 


80 



Chapter Five: Parameter Estimation mi Prediction of Uncertainly . 


’ z l(*l) 




■*l“ 

z 2 (t 2 ) 

— 


+ 

0 2 

_ z n{^n)_ 




_V 


or (5.1.1-4) 

z = f(t;cc) + 6 

Since the measurement error, 6, is assumed normally distributed, the 
measurement vector is also normally distributed with the following statistics: 

p = E{z} = E{/(t;a)} + E{fl} = /((;«) (5.1.1-5) 

e{(z - p){z - ft ) T } =E{ee T } =e ( 5 . 1 .i-6) 


5.1.2 Maximum Likelihood Estimation 

The maximum likelihood estimate selects the parameters a such that 
the probability of z, the likelihood function, is maximized. This method 
estimates parameter values so that the "truth" model observations are the 
most likely ones to have occurred. Since z is a normal random variable, its 
probability density is the normal joint probability density: 

r(z-,‘ a) = ‘ N - ex p[-M z ~-7(f;«)) T ©-‘(5 - /(»;*))] 

The maximum likelihood estimate maximizes the likelihood function 
(5.1.2-1) and minimizes the negative log-likelihood function below: 

£(z;a) = -ln[p(z;a)] (5.1.2-2) 

f(z ;a) = [f ln(z*) + iln(|0|)] + i ( z " (z “ /«;«)) 

(5.1.2-3) 

The maximum likelihood estimate Sc therefore satisfies the equation 
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(5.1.2-4) 


The maximum likelihood estimate h is obtained starting with an 
initial guess of the parameters, a 0 . Based upon the theoretical observations 
generated with a Q , adjustments, Aa, to the guessed parameters are sought 
such that 


A a = a-a 0 (5.1.2-5) 

Expanding Equation (5.1.2-4) in a Taylor series expansion about a 0 yields 


- -\\T 


0 = 


<?C(z;q Q 

da 


) 


-\\T 


a=a n 


d f d£fra) 
<?av da 


) 


Aa + O^Aa 2 ) (5.1.2-6) 


a=a n 


The parameter adjustments, Aa, are determined by taking the partial 
derivatives of the negative log-likelihood function from Equation (5.1.2-3) 
(assuming that the term in brackets does not depend on the parameters a) 
and by neglecting the higher order terms in the Taylor series expansion. The 
second derivative of the negative log-likelihood function with respect to a is 
replaced with its expected value, which can be theoretically expressed as the 
expected value of the dyadic product of first partial derivatives [42]: 


da 


r dj(z;a) ^ 

da 




= E 


<?C(z;or) 

da J da 


(5.1.2-7) 


Assuming that 0 does not depend on a, Equations (5.1.2-3) and (5.1.2-7) imply 
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(5.1.2-8) 


Replacing the Hessian of second partial derivatives in Equation (5.1.2-6) by its 
expected value (the Fisher information approximation) leads to the linear 
matrix equation, called the normal equations, for the adjustments to the 
parameters. These equations are formed with the partial derivatives of the 
theoretical observations with respect to the estimated parameters, the 
observation residuals, and the measurement error statistics. 
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e-'(z-f(z;a 0 )) 

(5.1.2-9) 


5.1.3 Weighted Least Squares Estimation 

Least squares estimation selects the parameters 5 so as to minimize the 
sum of the squares of the deviations between the "truth" model and 
theoretical measurements. For N observations, the least squares estimate 
seeks to minimize the quantity 

Q=|>,-/<(»,; a)) 2 (5.13-1) 

i =1 


If the measurements are of varying quantities and units, and some 
measurements are more reliable than others, the weighted cost function Q is 
used. 


Q _ V 1 ( 2 « fi(*i' a )) 

^ ~ 2-1 U7 

i=1 r Y i 


(5.1.3-2) 


where Wj is a sequence of weighting values [21]. A sequence of observations 
may involve a wide variety of physical quantities, i.e. distances, angles, 
temperatures, frequencies. Typically the error weightings non- 
dimensionalize these disparate measurements in the cost function. 

The least squares or weighted least squares estimate is found by setting 
the derivative of the cost function with respect to the estimation parameters 
a equal to zero. 
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5.1.4 Choice of Estimation Method 


The weighted least squares and maximum likelihood methods provide 
the same estimates if the measurement residuals, 6, are uncorrelated. 

E{ 0, Gj } = 0 V ijtj (5.1.4-1) 


By assuming that this measurement noise is uncorrelated, the covariance 
matrix, 0 defined in Equation (5.1.1-6) becomes the diagonal matrix 


of 0 
0 g\ 


0 

0 



(5.1.4-2) 


where is the variance and a n is the standard deviation of the n th 
measurement. Inserting this diagonal matrix into Equation (5.1.2-1) and 
(5.1.2-3) for the likelihood and negative log-likelihood functions results in the 
following: 


p(z;a) = 
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(2k) 2 g 1 -g 


■exp 
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1 V 1 ( z i 


-il 


i=l 


Gf 


(5.1.4-3) 


C ( 2 ; a ) = [f ln(2«) + ln(<Ti • ■ • a N )]+^ ^ 


(5.1.4-4) 


Both the maximum likelihood and weighted least squares estimates 
are determined by setting the partial derivative of a cost function with respect 
to a equal to zero. Since the maximum likelihood estimate assumes that the 
constant part in brackets [ ] from (5.1.4-4) does not depend upon the 
parameters a, the two methods (compare Equations (5.1.3-2) and (5.1.4-4)) 
provide identical estimates when uncorrelated measurement error variances, 
cl, are used to weight the measurement residuals. Since PEP's parameter 
estimation routine assumes uncorrelated measurement errors, its estimation 
technique is referred to as least squares maximum likelihood estimation. 
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5.1.5 Incorporation of A Priori Information 

Frequently, existing models provide a priori information about some 
of the estimation parameters. The estimation process can be shortened or 
simplified by including this a priori information in the normal equations. 

Suppose a priori information has estimated the value of the parameter 
a j as a, with an uncertainty, or standard deviation, of d,. These a priori 
parameter values are then grouped into the m-dimensional vector cc, where 
m is the total number of parameters estimated. A zero value is assigned to 
any parameters when no a priori knowledge is available. The variances of 
the a priori estimates are collected into the m x m diagonal matrix X. For any 
parameters without a priori standard deviation information, their diagonal 
element is infinity, although in practice some reasonably large value will 
suffice. A full covariance matrix can be used for Z if correlations are 
available for the a priori parameter estimates a. 

With a priori information, instead of minimizing Equation (5.1.2-3), 
the following functional is minimized [11]. 

Q =(z - 7(t;fi))V (z - /((; &))+(a - a) T ±-' (a - a) (5.1.5-1) 

For parameters where no a priori standard deviation information is available, 
the cost function remains the same as in (5.1.2-3). 


Setting the partial derivative of this cost function with respect to the 
parameters a equal to zero and making the same assumptions as in Section 
5 .1.2, yields the following linear matrix equations. 


A = 




da 
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0- 1 


df(z;a 0 ) | e.-, 
da 


(5.1.5-2) 


b = 


f df(z;a 0 ) ^ 
da 


0-‘ (z - 7(2; & 0 )) - XT 1 (a„ - fi) (5.15-3) 


such that 

A a = A’ 1 b (5.1.5-4) 
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From studies of the earth's gravitational field, William Kaula proposed 
an empirical relationship for the variance of the gravitational harmonic 
coefficients [24]. Estimates of the lunar gravitational harmonic coefficients led 
him to believe that the same relationship was valid. For the lunar 
coefficients, Kaula adjusted the constant coefficient to account for the 
gravitational field strength difference between the moon and earth [37]. From 
their studies of the lunar gravitational field. Bills and Ferrari chose a slightly 
different semi-empirical formula for the covariance; one which had been 
proposed by Vening-Meinesz. After having difficulty converging upon 
gravitational harmonic coefficients. Bills and Ferrari included the Vening- 
Meinesz a priori covariance information and developed a 16 x 16 lunar 
spherical harmonic potential model with coefficient uncertainties [12]. 
Research scientists at JPL continue to use this method to develop 50 X 50, and 
even 75 X 75 spherical potential models of the moon [27]. 

The method of incorporating a priori information into the normal 
equations was not used for these simulations. This thesis evaluates 
gravitational sensing schemes which may be employed to develop gravity 
field models. As such, an arbitrary lunar gravitational "truth" model was 
developed and there was no a priori information available regarding this 
"truth" model. When a gravitational sensing method is selected and its 
mission flown, a priori lunar gravitational field information can and should 
be used in the development of new gravitational field models. 


5.1.6 Solution to the Normal Equations 

Since the normal equations [(5.1.2-9) or (5.1.5-4)] are linear, the 
adjustments to the parameters can be solved by various numerical 
techniques. PEP uses the Gauss-Jordan method to simultaneously solve the 
normal equations and invert the coefficient matrix. This method uses 
diagonal pivots without interchanging rows or columns, so that only the 
lower diagonal half of the symmetric coefficient matrix A is stored in 
memory. 

Numerical problems can arise if the coefficient matrix is ill- 
conditioned. This could happen if the vector of parameter adjustments A a 
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consisted of widely different sized quantities or quantities with widely 
different units. To prevent these numerical difficulties, PEP scales the 
normal equations prior to solving them. This process scales each row and 
column of the symmetric coefficient matrix by the square root of the diagonal 
element. To preserve the equality, each row of the right hand side vector is 
also scaled by the same factor. After the adjustments to the parameters are 
solved, PEP unscales its rows to provide the proper units and values to the 
parameter adjustments [8]. 

This scaling process can be avoided by using the square root of the 
normal equations [11]. This method takes advantage of the properties of the 
symmetric coefficient matrix and, by using the square root of the normal 
equations, lessens the effect of disparate units or scale factors. The JPL orbit 
fitting software uses this method [27]. 

Since the equation for A a was obtained by neglecting second and 
higher order terms in a Taylor series expansion, this adjustment will not 
yield St exactly. An iterative technique must be used to approach a maximum 
likelihood estimate for the parameters. Once the parameter adjustments are 
determined, the initial parameter guesses are adjusted. The equations of 
satellite motion are then re-integrated and the theoretical observations, 
partial derivatives, and residuals are recalculated. The normal equations are 
then reformed and a second set of parameter adjustments are determined. 
These iterations continue until the process converges. 


5.1.7 Other Estimation Techniques 

For some applications, least squares maximum likelihood may not be 
the optimum method for estimating the gravitational coefficients. If noise 
were included in the satellite dynamics, Kalman filtering techniques would 
be more appropriate. A linearized Kalman filter about a nominal satellite 
trajectory or an extended Kalman filter for which the reference trajectory is 
updated for each observation could be used [31, Chapter 9]. The Kalman filter 
would then estimate the gravitational parameters as well as the satellite's 
state by augmenting the state vector to include them as non-dynamic state 
variables. 
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Since Kalman filtering techniques become extremely cumbersome as 
higher degree models are estimated, a more efficient technique may be 
necessary. Maximum likelihood system identification combines the two 
approaches by performing a Kalman filter on the satellite state, X*, and a 
maximum likelihood estimate on the satellite's orbital initial conditions and 
the gravitational harmonic coefficients [31, Chapter 10]. In this method, the 
satellite's motion is propagated with a Kalman filter and then a maximum 
likelihood adjustment is made to the initial conditions and parameters. The 
process is then repeated until convergence. 

Unmodeled forces such as radiation pressure, fuel tank leakage, and 
higher order gravitational harmonic effects will cause noise in the satellite 
dynamics. Because of this noise, one of these techniques should be used. 
Least squares maximum likelihood estimation, however, was sufficient for 
this thesis' evaluation of sensing methods. 


5.2 Statistical Prediction of Uncertainty 

Solving the normal equations produces information about the 
uncertainty of the estimated parameters St as well as the values of the 
parameters themselves. The parameter estimation covariance matrix 
provides a measure of the uncertainty of the estimates. This additional 
information can be used to evaluate the resulting estimated model of the 
observed behavior. 

The coefficient matrix A from (5.1.2-9) and (5.1.5-4) is the Fisher 
information matrix and by the Cramer-Rao lower bound, the covariance of 
any unbiased estimate is greater than or equal to its element in the inverse of 
the Fisher information matrix [42, 46]. Additionally, maximum likelihood 
estimates have the desirable property that they are asymptotically unbiased 
with the Cramer-Rao lower bound obtained as the number of observations 
approaches infinity [42]. This bound only applies if the nonlinear system has 
been modeled correctly in (5.1.1-4). 
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If the assumed standard deviation of the measurement errors, a„, are 
correct and the modeling is correct for both the equations of motion and the 
observations, the root mean square (rms) of the post-fit observation residuals 
divided by the assumed measurement errors would be approximately one: 

rms = ^N^(z-/(*;a)) 0"‘(z-/(f;a)) (5.2-1) 

Typically the number of observations is much greater than the number 
of parameters, so the 1 / (N-m) term is replaced by 1 /N. When the above rms 
differs greatly from one, the parameter covariance matrix produced by the 
inverse of the coefficient matrix from the normal equations (A' 1 ), by the 
Cramer-Rao lower bound, should be multiplied by this rms to obtain a truer 
estimate of the uncertainty in the parameter estimates. This adjustment 
accounts for incorrect values of the measurement standard deviations, c n in 
0 , but does not account for any modeling error effects on the uncertainty. 

The parameter covariance matrix, X, is then obtained as 

I = rms x A' 1 (5.2-2) 


The standard deviation of parameter estimate a, is 

cr,. 

and the correlation between parameter estimates a, and ccj is 


a,<Jj 


(5.2-3) 


(5.2-4) 


PEP saves the matrix A' 1 and the rms of the observation residuals 
divided by the assumed measurement errors resulting from an estimation 
run so they can be used to predict the uncertainty of an orbit propagated with 
the estimated parameters a. 

Let X* be the state of a lunar satellite at time to assumed to be known 
perfectly, perhaps by determination from navigation satellite observations. 
These satellite initial conditions are numerically integrated in time using the 
estimated gravitational model. Partial derivatives of the satellite's motion 
with respect to the parameters a are also numerically integrated. The state 
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covariance matrix, S, for the uncertainty in the satellite state, can then be 
calculated at any time t as 



dXjo 

da 


Z(t) 


f dX'm 
K da 


\ T 


(.5.2-3) 


The standard deviation of the uncertainty in any one coordinate 
direction is the square root of that diagonal element. The root sum square 
(rss) of the position and velocity uncertainties are then 

rss(-R) = ^/ijj + i 2 2 + “33 (5.2-4) 

rss( V ) = -\]S 44 +H 55 +H 66 (5.2-5) 


Since the state covariance matrix 5 is based upon the augmented state 
vector (3.3-1), it is partitioned as follows: 



E m (t) 



(t) 

in 


(5.2-6) 


Often, the cross correlations between position and velocity E RV are neglected 
and the covariance matrix for position E RR and velocity E vv are used 
separately to evaluate the uncertainty of orbit prediction. The state 
uncertainty analyses performed for the estimated lunar gravity fields in 
Chapter Seven use these individual state covariance matrices rather than the 
augmented state covariance matrix. 

The inertial (x,y,z) coordinate uncertainties thus obtained can also be 
converted to a local vertical, local horizontal coordinate frame to provide 
navigation uncertainty information. The local vertical, local horizontal 
frame is defined by the vertical (VT), down range (DR), and cross track (CT) 
directions. Their translations from inertial coordinates are obtained by 


u VT (t) = unit[X(t)] 

(5.2-7) 

u CT (o = umf|X(o x Xu)] 

(5.2-8) 

Udr (0 = HyT (0 ^ M qy (0 

(5.2-9) 
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The transformation matrix from local vertical, local horizontal to 
inertial coordinates is then 

Rjl (0 = i UQj(t) i ^dr (5.2-10) 

and this transformation can be used to obtain the covariance matrix for the 
uncertainty in the satellite state in the LVLH frame. 

Sf<0 = Ri«»2“(f)R IL <„ (5.2-11) 

3r«> = KhoZVit) Ril(') (5.2-12) 

The standard deviation of the uncertainty in any one LVLH coordinate 
direction is the square root of that diagonal element. With these formulas, 
the uncertainty of position and velocity in the vertical, cross track, and range 
directions can be calculated at each time t due to uncertainties in the 
estimated parameter a. These routines were coded in analytical software 
separate from PEP and each estimated gravity model was evaluated in this 
manner. This analysis was performed for both the PEP supplied covariance 
matrix (A -1 ) and the covariance adjusted by the rms of the observation 
residuals divided by the assumed measurement errors. These analyses 
simulated the propagation of uncertainties onto the far side of the moon after 
obtaining a navigation fix on the near side. 

The previous statistical uncertainty prediction was performed using 
the assumption that the lunar satellite initial state, X 0 , was known perfectly. 
If this is not the case, the same process for calculating the uncertainty of orbit 
prediction would still apply. In general, the vector of estimated parameters oc 
would be augmented to include estimates of the satellite initial conditions 
and the process repeated. Since this thesis is concerned with the uncertainty 
due to a mismodeled lunar gravitational field, perfect knowledge of the 
satellite's initial conditions was assumed. 
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5.3 Prediction Uncertainty Due to Mismodeling 

The statistical uncertainty analysis presented in the previous section is 
based upon the assumption that the estimated gravitational model correctly 
models the lunar gravity field. Due to the inclusion of mascons in the 
"truth" model, the finite spherical harmonic estimation model could not 
correctly model the observed behavior. In most real world cases there is no 
way to account for or analyze the errors between modeled and true behavior 
because the modeling errors cannot be separated from the other errors. For 
this thesis, however, a direct comparison between the true and predicted 
behavior can be made, since both the "truth" and fit models are available. 

Given satellite initial conditions, X*, orbits for the truth and fit models 
can be numerically integrated. At each point in time, the position and 
velocity deviations between the two models can be calculated. These 
modeling errors in position and velocity can then be transformed from 
inertial coordinates to the local vertical, local horizontal frame using the 
transformations (5.2-7), (5.2-8), and (5.2-9). The root sum square of the 
position and velocity errors due to gravitational field mismodeling can also 
be calculated. The impact mismodeling has upon navigation uncertainty 
prediction is revealed by comparing the modeling errors with the statistical 
uncertainty predictions. 

Once again, using software separate from PEP, this analysis was 
performed for orbits propagated from the lunar near side to the far side for 
one orbital period. Additionally, the uncertainty due to mismodeling was 
analyzed by comparing the "true" path of a lunar landing when the 
maneuver was calculated based upon simulations using the estimated model 
and then executed in the "true" lunar gravitational field. These landing 
errors provide a feel for the impact lunar gravitational field mismodeling 
will have upon the execution of future space missions. 
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Orbit Fitting 
and 

Gravity Estimation 


6.1 Methodology 

Several different lunar gravitational sensing methods were evaluated 
using the PEP orbit integration, observation generation, and parameter 
estimation capabilities described in Chapters Three, Four and Five. For each 
sensing scheme, the following procedures, depicted in Figure 6.1-1, were used 
to generate "truth" model observations and estimate gravitational field 
coefficients. 

First, each sensing scheme was analyzed to determine the nature of the 
satellite orbits employed, the type of observations generated, and the 
observation accuracy and frequency. Using the "truth" model and these 
descriptions of the observation method, PEP simulated the mission's orbits 
and observations, including lunar occultations. Auxiliary software routines 
converted the integration output to selenographic osculating orbital elements 
(a,e,i,£2,co,M) as well as selenographic spherical coordinates (r,0,<t>) versus time. 
This data was plotted to analyze the orbit's stability over the integration 
period (14 to 28 days). For stable orbits, the "truth" model observation file was 
then used to estimate the fit model lunar spherical harmonic coefficients and 
satellite initial conditions. 

For the estimation process, the gravitational harmonic coefficients and 
satellite initial conditions were perturbed from the "true" values used to 
generate the observations. The gravity field was altered to reflect the fit 
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model's degree and order ("truth" model usually included mascons). From 
the initial perturbed parameters, a 0 , the satellite orbits were propagated, 
theoretical values of the observations were calculated, and the difference 
between the "truth" model observations and these theoretical observation 
values (observation residuals) were calculated. Partial derivatives of the 
motion and partial derivatives of the theoretical observations were calculated 
with respect to all parameters to be estimated. From the observation residuals 
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Figure 6.1-1: Flow Diagram for Gravity Sensing Mission Simulations 

and partial derivatives, the parameter estimation module formed and solved 
the normal equations to calculate the adjustments to the parameters, A a. 
The size of the parameter adjustments and their uncertainties determined if 
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the process had converged. If the process had not converged, the new 
parameters were used to propagate the next iteration’s orbit and the process 
was repeated. The maximum likelihood parameter estimates and the 
covariance matrix (A -1 ) were written to output files when the process finally 
converged. These files were then used to analyze the estimated gravitational 
model and its ability to predict future lunar missions. 

Although this process is straight-forward, it is a very time and 
computer memory intensive process. A typical 14 day orbit propagation for 
an 8 X 8 spherical harmonic gravity field with partials calculated for 79 
parameters typically took over eight hours on a Sun Sparcstation IPC. Also, 
due to the high correlations between gravitational parameters, it was not 
unusual for a run with different "truth" and fit models to require at least 
fifteen iterations to converge upon a solution. The convergence criterion, 
perhaps unnecessarily stringent, required the adjustments to be less than 1% 
of the parameter uncertainties. 


6.2 Fifth Degree Harmonic Truth and Fit Model Test Case 

Since PEP was modified to include lunar mascons and expand the 
degree and order spherical harmonic expansions allowable, several test cases 
were rim to ensure that PEP-D's modules were operating correctly. The first 
set of tests involved the attempt to estimate a 5 X 5 spherical harmonic 
expansion fit model for observations generated with a 5 X 5 spherical 
harmonic expansion truth model. Since the degree and order of the 
estimated fit and truth models were the same, this test was expected to be a 
good warm-up exercise for future runs. 

For these tests, a single lunar polar 200 km altitude near-circular 
satellite orbit was numerically integrated for twenty-eight days. The initial 
selenographic orbital elements are listed in Table 6.2-1, where elements in 
parentheses are in PEP's internal units (angles referred to the mean equinox 
and equator of the earth of 1950.0). A zero degree initial mean anomaly was 
selected first, placing the satellite directly over the lunar north pole. Due to a 
singularity in PEP's algorithm for the central body harmonic effects, it could 
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not propagate this orbit. Because of perturbation forces, it would be highly 
unlikely that this condition would be repeated during the middle of an orbit 
propagation The mean anomaly was adjusted one degree and the process 
restarted with no further difficulties. 

Table 6.2-1: Satellite Initial Conditions for 5 X 5 Spherical Harmonic Test Case 


ao 

1938 km (1.295472995 x 10 5 AU) 

eo 

0.05 

io 

90° (103.1048493849350°) 

Do 

90° (304.1996805394721°) 

0)0 

90° (72.2394269798987°) 

Mo 

1 ° 


The starting time for the orbit integration was 16 May 1968 or Julian 
Date 2,440,001.5 0 hour Coordinate Time. This initial epoch was selected over 
an epoch in 1996 or 1997, when a lunar gravitational sensing mission might 
occur, because it is near the beginning of the n-body file supplied by the SAO. 
This file is read and interpolated from during the integration and observation 
modules, so computer time was saved by selecting an initial epoch near the 
beginning of the file. Other than saving computer time, the choice of initial 
epoch should not effect the simulations of this thesis. 

The lunar satellite's orbit was propagated using the first five degree and 
order coefficients from the 1980 Bills and Ferrari 16 X 16 lunar spherical 
harmonic model [12]. These coefficients are given in Table 6.2-2. The central 
body attraction of the moon and the perturbing attractions of the earth and 
sun were included in the Adams-Moulton numerical integration with a step 
size of 2 -14 day. This allowed approximately 1,451 steps per 2.125504 hour 
orbit. This was an extremely small step size to characterize fifth degree and 
order harmonics and ensure numerical stability in the integration of motion. 

The observation module then created Doppler observations of the 
satellite from NASA DSN's Goldstone site every minute. Four different 
observation scenarios were generated for this test case. Two involved range 
and Doppler observations and two involved Doppler observations only. In 
two cases observations were occulted by the moon and in the other two cases 
the observations continued through a fictitious "transparent" moon. 
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Table 6.2-2: 5X5 Spherical Harmonic Truth Model Coefficients 


Harmonic 

xlO- 6 

Harmonic 

xlO- 6 

h 

202.431 



C 2 i 

-0.07 

S 21 

-0.00 

C 2 2 

34.49 

S 22 

0.03 

J3 

8.8897 



C 3 1 

21.96 

S 31 

6.63 

C 32 

14.14 

S 32 

4.76 

C 33 

15.87 

S 33 

-2.45 

J4 

-11.73 



C<n 

-4.82 

£>41 

1.91 

C 42 

-8.13 

S 42 

-6.76 

C 43 

0.48 

S 43 

-14.43 

C 44 

-3.50 

S 44 

-0.55 

J5 

2.388 



C 5 1 

-9.66 

S 51 

-1.53 

C 52 

3.71 

S 52 

-2.35 

C 53 

-0.39 

S 53 

4.91 

C 54 

0.56 

S 54 

-6.58 

C 55 

-6.69 

S 55 

11.60 


The initial test cases recreated the measurements of the Lunar Orbiter 
missions (Section 1.3). The orbits were selected to simulate orbits proposed 
for the Lunar Observer or Lunar Scout missions [16, 25, 39]. The observation 
scheme, duplicating the processing of Lunar Orbiter data, used the DSN S- 
band frequency (2.115 GHz) for collecting Doppler range rate observations for 
sixty-second intervals. No transponder frequency translation was simulated. 
These observations had a quoted accuracy of 1 mm/sec [29], which is the 
accuracy achievable for integer Doppler counts over the sixty-second interval. 
Range observations were also generated in the "truth" model observation 
file. When fitting to only Doppler observations, the range observations were 
given negative error weights so that they would be ignored in the estimation 
process. The two way range measurement was included on the same signal 
used for the Doppler observable and had a three meter accuracy. This 
accuracy is based upon the ability to detect a 20 nanosecond time delay for a 
two-way phase shift keyed code. 
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Over the twenty-eight day observation period, the satellite would orbit 
over the entire lunar surface as the moon rotated under its orbital plane. Due 
to the extremely conservative integration step size, the twenty-eight day 
viewing period, and the calculation of partials, these initial test cases quickly 
filled up the available disk space, resulting in a system crash. Reevaluating 
the simulation requirements, the observation period was reduced to fourteen 
days, which still provides complete lunar coverage in either an ascending or 
descending pass. 

After the “true" observations were generated, the satellite initial 
conditions and thirty-two lunar harmonic coefficients were perturbed to 
provide the first guess for the estimation process. Each gravitational 
harmonic coefficient's absolute value was increased by 1.0 x 10' 7 . The orbital 
initial conditions were altered by the values given in Table 6.2-3. The satellite 
orbit was then numerically integrated along with partial derivatives. 

Table 6.2-3: Perturbed Satellite Initial Conditions for Estimation Model 


8a 

149.6 km (1.0 x 10' 11 AU) 

8e 

1.0 xlO- 5 

8i 

-1.0 xlO- 4 0 

5 Q 

-1.0 x 10-4 ° 

So) 

-1.0 x 10- 4 0 

8Mq 

1.0 x 10- 4 0 


These tests were run using saved partial derivatives of motion. Using 
this feature, partial derivatives were calculated during the first iteration step 
and these saved partials of motion were used to determine the partial 
derivatives of observations on subsequent iterations. Computation time was 
saved by not recalculating the partials on each iteration. This approach 
worked well for test cases with the same "truth" and fit models, but was not 
effective on other simulations. In general, this feature is very effective when 
the process is close to convergence. 

The fitting process revealed that all of the harmonic terms included in 
the equations of motion have to be included in the partial derivative 
calculations. By default, PEP was only calculating the second degree harmonic 
effects upon the partials. This provides sufficient accuracy for earth satellites 
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where the ]i effect dominates motion, but it was not sufficient for lunar 
satellites. Additionally, the first scenario's slow convergence raised concerns 
about the choice of initial osculating orbital elements. For a circular orbit, the 
argument of perilune (to) is undefined. For the near-circular orbit, the 
perilune position was very difficult to determine. The argument of perilune 
and mean anomaly were difficult to estimate because of their high 
correlation. PEP software does not currently use equinoctal orbital elements 
[e sin(Q + co), e cos(D + co)], but it provides the option of replacing the Q, to, and 
M orbital elements with the sum of the angle elements £2, Q + co, Q + co + M. 
After the initial orbital elements were converted to this form, their 
correlations were reduced and the parameter estimation routine had no 
difficulty estimating these parameters. 

After overcoming the previously mentioned difficulties, the 
simulations converged in three iterations for all four scenarios. After these 
test cases, a fifth scenario was run in which radar biases for the time delay and 
Doppler shift observations were estimated along with the gravitational 
parameters and satellite initial conditions. This scenario evaluated the 
estimation procedure's ability to handle measurement bias estimates and also 
converged in three iterations. 


6.2.1 Scenario One: No Occultation. Range and Range Rate 

The first scenario calculated a spherical harmonic fit model based on 
20,156 time delay (range) and 20,156 Doppler shift (range rate) observations 
over the fourteen day period. Observations were processed every sixty 
seconds for the entire period since there were no horizon or occultation 
restrictions. The thirty-two gravitational harmonic coefficients and six initial 
conditions were all estimated within approximately ten digits of accuracy for 
each parameter (six digits for angular initial conditions). The observation 
residuals divided by the assumed measurement errors, referred to as the non- 
dimensional "fit residuals", had a root mean square (rms) value of 7.45263 x 
10* 4 , essentially zero. The theoretical rms for these cases is zero, as opposed to 
the value of one mentioned in Section 5.2 because the "true" observations 
were exact (without measurement noise) and the truth and fit models are of 
the same degree and order. The run also calculated correlations as high as 
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-0.998 between the C 33 and C 53 coefficients and S 32 and S 52 coefficients. 
Fifteen of the 703 correlations were greater in magnitude than 0.95. Initially, 
these high parameter correlations seemed as if they were due to a resonance 
in a particular harmonic frequency, but later fits to higher degree models also 
demonstrated high correlations for the higher degree terms. 


6.2.2 Scenario Two: No Occultation, Range Rate 

This scenario excluded the 20,156 time delay (range) observations from 
the previous run's fourteen day period in its estimation of the spherical 
harmonic fit model. Once again the thirty-eight parameters were all 
estimated with approximately ten digits of accuracy and the fit residuals were 
essentially zero, with an rms value of 8.83418 x 10‘ 4 . The highest parameter 
correlation was again -0.998 between C 33 and C 53 . The S 33 S 53 and S 32 S 52 
correlations were the second highest at -0.997. The S 33 and S 53 correlation 
was up from -0.996 in the previous scenario. Once again fifteen of the 703 
correlations were greater in magnitude than 0.95. 


6.2.3 Scenario Three: Occultation, Range and Range Rate 

After the previous scenarios were completed, the "truth" model 
observations were recreated with lunar occultations. This third scenario then 
calculated its fit model based upon 15,672 time delay (range) and 15,672 
Doppler shift (range rate) observations over the fourteen day period. Lunar 
occultations eliminated 3 days, 2 hours and 44 minutes of observations over 
the fourteen day period. Again the thirty-eight parameters were estimated to 
within ten digits accuracy of their "true" values. With the reduced 
observability, fit residuals increased, having an rms value of 2.54901 x 10 -3 . 
Although these residuals are larger, the fit model is still essentially exact. The 
parameter correlations for this scenario essentially repeat those of the first 
scenario with the C 33 C 53 and S 32 S 52 coefficients having the highest 
correlation (-0.998). As in both previous cases, fifteen of the 703 correlations 
were greater than 0.95. 
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6.2.4 Scenario Four: Occultation, Range Rate 

This scenario excluded the 15,672 time delay (range) observations from 
the previous run in its estimation of the gravitational coefficients and orbital 
initial conditions. The scenario's fit model agreed with the truth model to 
about nine places. The fit residuals for this scenario were up from the 
previous scenario with an rms value of 6.86066 x 10' 3 . The highest parameter 
correlation was the C 33 C 53 correlation (-0.999). The S 32 S 52 and S 33 S 53 
correlations were once again close behind (-0.998). As before, fifteen of the 703 
correlations were greater than 0.95. 


6.2.5 Scenario Five: Occultation, Range and Range Rate with Biases 

This scenario used the 15,672 time delay (range) and 15,672 Doppler 
shift (range rate) observations from the third scenario and attempted to 
estimate the thirty-two gravitational harmonic coefficients and six orbital 
initial conditions as well as measurement biases for the time delay and 
Doppler shift observations. This scenario estimated the gravitational 
coefficients and initial conditions with about nine places of accuracy. 
Additionally it estimated a 1.43121 x 10 ' 9 second time delay bias and a -4.37639 
x 10 ' 9 Hz Doppler shift bias. No biases were included in the observations and 
this run essentially estimated zero biases. The fit residuals were lower than 
the fourth scenario's with an rms value of 6.10319 x 10‘ 3 , also essentially zero. 
The C 33 and C 53 coefficients once again had a -0.999 correlation with the 
S 32 S 52 and S 33 S 53 parameters yielding the second highest correlations (-0.998). 
Fifteen of the 780 correlations were greater than 0.95. 


6.3 Earth-Based Doppler Observable Mascon Test Cases 

The second set of tests included two small mascons in the truth model 
and increased the degree of the spherical harmonic fit model. This test case 
was run to determine whether PEP-D could estimate a spherical harmonic 
model based on observations with a different truth model and to determine 
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whether the sensing geometry affected the estimate's ability to model the 
lunar gravitational field. 

For the first scenario, the previous 5X5 spherical harmonic expansion 
(Table 6.2-2) was augmented with two mascons of equal strength placed at the 
limbs of the moon on the equator. For the second scenario, the 5 X 5 model 
was augmented with mascons of equal strength placed on the front and back 
faces of the moon, still at the equator. Mascons were included in pairs to 
preserve the center of mass of the moon and avoid solving for first degree 
harmonic coefficients. The strongest surface disk estimated in the Wong 
model [47] was selected for the strength of these mascons whose strengths and 
locations are listed in Table 6.3-1. 


Table 6.3-1: Mascon Test Case - Mascon Placement 



Mascon on Limb 

Mascon on Face 

1 

2 

1 

2 

m, 

7.212 x 10- 6 

7.212 x 10- 6 

7.212 x 10- 6 

7.212 x lO* 6 

r i 

1680 km 

1680 km 

1680 km 

1680 km 

mm 


1 



1 b 






The same lunar polar 200 km altitude near-circular satellite orbit was 
numerically integrated for fourteen days over these two lunar gravitational 
fields. For both cases, the orbit integration was once again started on 16 May 
1968 and used the same satellite initial conditions as in the previous tests 
(Table 6.2-1). The lunar polar satellite's motion was detected by the same 
earth-based Doppler shifted observables used for the previous tests. 
Observations were again interrupted by lunar occultations. To reduce the size 
of the computer files generated in the simulation, the Doppler count interval 
was changed from sixty seconds to one hundred and twenty seconds. 

These tests were run to determine whether earth observations are 
sufficient to estimate the lunar gravitational field despite the mascons' 
location. Mascons cause radial disturbing accelerations as satellites pass over 
them and Doppler observations sense relative velocity along the line-of-sight. 
Doppler observations should therefore have an easier time sensing mascon 
disturbances when the observing site, observed body, and mascon are all 
aligned. These two lunar mascon scenarios present two geometrical extremes 
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for Doppler sensing. In the first case, the disturbing acceleration due to 
passage over a mascon on the limb is orthogonal to the Doppler line-of-sight. 
In the second case, the disturbing acceleration and Doppler line-of-sight are 
aligned as the orbiting satellite passes over the near face's mascon. 

These test cases were initially run using the saved partials feature of 
PEP mentioned in Section 6.2. For the first scenario full partials were 
recalculated on the fifth iteration and the parameters converged on the ninth 
iteration. The parameters from this solution were then run with full partials 
and they required six more iterations to converge. For the second scenario, 
the first five iterations were run using the saved partials feature since the 
"semi-convergence" of the previous case had not been discovered. After the 
fifth iteration, the saved partials feature was abandoned and this scenario also 
required fifteen total iterations to converge. 

The mascon test cases, requiring fifteen iterations to converge, 
demonstrated the difficulty of estimating gravitational parameters when the 
"truth" and fit models differed. Additionally, these scenarios revealed the 
limitations of PEP's saved partials feature. Since it was not clear that saved 
partials were helping the estimation process, the method was abandoned. 
This significantly increased the amount of computer time and the size of the 
memory files required for the simulation process. 


6.3.1 Scenario One: Limb Mascons 

As mentioned previously, this scenario converged in fifteen iterations 
to a degree and order eight spherical harmonic model. Appendix E's Table E-l 
lists the estimated harmonic coefficients for this fit model. Additionally this 
run estimated the satellite osculating orbital initial conditions. Table 6.3.1-1 
compares the true initial conditions with their estimated values. 

The fit residuals (observation residuals divided by assumed 
measurement noise) for this model have an rms of 7.36725. For the first time 
this non-dimensional statistic, which measures how well the estimated 
model agrees with the observed behavior, is greater than one. This shows 
that this spherical harmonic model is inefficient at modeling a gravitational 
field with mascon anomalies. An analysis of the parameter correlations 
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reveals that four pairs of correlations have cross correlations greater than 
-0.999: S66^86/ C66C86/ C55C75, and S55S75. The C54C74 and S54S74 

correlations also had an extremely high correlation of -0.989. Twenty-two of 
the 3,403 correlations were greater in magnitude than 0.90 with nine of these 
greater than 0.95. These correlations for the same 200 km lunar polar orbiter 
provided an indication that the high parameter correlations discovered in 
Section 6.2 were not due to the resonance of a particular harmonic frequency. 


Table 6 . 3 . 1 - 1 : Limb Mascon Orbital Initial Condition Estimates 



"True" 

Fit 

A 

ao 

1938.0 km 

1938.1106 km 

110.6 meters 

eo 

0.05 

0.050091 

0.000091 

io 

103.1048° 

103.1067° 

0.0019° 

Qo 

304.1997° 

304.1939° 

-0.0058° 

(D-ko)o 

16.4391° 

16.6468° 

0.2077° 

(£2-kd+M)o 

17.4391° 

17.6389° 

0.1998° 


6.3.2 Scenario Two: Face Mascons 

This scenario fit the observed behavior to an 8 X 8 spherical harmonic 
expansion in fifteen iterations. Table E-2 in Appendix E lists the estimated 
harmonic coefficients for this fit model. This run also estimated the satellite 
orbital initial conditions for the observed lunar satellite. .Table 6.3.2-1 
compares the true initial conditions with the estimated initial conditions. 

This case's non-dimensional fit residuals had a root mean square of 
11.8025, once again an order of magnitude greater than desired for an 
estimated model. The increase in the residuals' rms over the previous 
scenario is most likely due to the lack of observations for the far-side mascon. 
The Doppler observations in this scenario sense the radial accelerations due 
to the mascon on the near face but cannot account for the lunar far-side 
perturbations with the 8x8 spherical harmonic expansion. In the previous 
case, since the mascon disturbing accelerations were orthogonal to the 
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Doppler line-of-sight, the observations did not sense large disturbances and 
had an easier time fitting the small disturbances to the fit model. 

An analysis of the parameter correlations once again reveals that four 
pairs of correlations are greater than -0.998: S 66 S 86 / C 66 Cs 6 r C 55 C 75 , and 

S 55 S 75 . The C 54 C 74 and S 54 S 74 correlations are again very highly correlated 
at -0.989 and -0.988 respectively. Of the 3,403 correlations, nineteen were 
greater in magnitude than 0.90 and nine of these were greater than 0.95. 


Table 6.3.2-1: Face Mascon Orbital Initial Condition Estimates 



"True" 

Fit 

A 

ao 

1938.0 km 

1937.9781 km 

-21.9 meters 

eo 

0.05 

0.049956 

-0.000044 

io 

103.1047° 

103.1015° 

-0.0032° 

Qo 

304.1997° 

304.2022° 

0.0025° 

(£2+cd)o 

16.4391° 

16.3311° 

-0.1080° 

(D+gh-M)o 

17.4391° 

17.3624° 

-0.0767° 


6.4 Truth Model Development 

Since the true lunar gravitational field is not precisely known, a lunar 
gravitational "truth" model was developed for this thesis. This truth model 
was then used to evaluate various proposed lunar gravitational sensing 
schemes. Previous lunar missions have discovered mascons on the near side 
of the moon [35]. There is no conclusive evidence regarding mascons on the 
far side of the moon, since there are no observations of a satellite's motion 
over the far side. For this thesis' truth model, mascons were placed on both 
the near and far sides. These mascons augmented the previously used 5x5 
spherical harmonic expansion (Table 6.2-2). 

From their surface layer representation of the lunar gravitational field, 
Wong et. al. identified major lunar mass anomalies with their selenographic 
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features [47 Table 2]. Wong's group modeled the lunar gravitational field 
with 50 km radii surface disks distributed around the lunar surface (Section 
2.3.2). Each surface feature was therefore represented by several disks 
depending upon the size of the feature. The seven most influential mass 
anomalies, requiring forty-nine surface disks, were selected for the truth 
model's near-side mascons. Since PEP-D does not model mascons as surface 
disks, these forty-nine disks were converted to point masses placed 58 km 
below the lunar surface. The selenographic positions and strengths of the 600 
surface disks [47 Figure 4 and Table 6] were then correlated with topographic 
and gravitational maps depicting the major surface features [26] to determine 
the forty-nine point masses required for the "truth" model. 

Since all of the previously identified mascons are on the near side, 
there was no scientific basis for establishing lunar far-side mascons. For this 
thesis, far-side mascons were developed to meet two requirements. First, they 
should be difficult to detect from earth-based observations, and secondly they 
should not alter the lunar center of mass. Since a mascon placed on the back 
face of the moon would be the most difficult to sense, the first far-side mascon 
was centered at 180° longitude and 0° latitude. Masses comparable to the 
near-side mass anomalies were selected for these point masses. Since lines of 
longitude are more widely spaced at the equator and nineteen point masses 
were used, this mascon is the largest and strongest of any mass anomaly in 
the truth model. 

The translation of the lunar center of mass determined the placement 
of the second far-side mascon. This balancing mascon was composed of seven 
point masses centered on 214.5° longitude and -48.0° latitude. After these two 
far-side mascons were added to the near-side mascons, the center of mass was 
still askew, so an additional point mass added to the model. The entire 
gravitational truth model is contained in Appendix D. 

The lunar polar 200 km altitude orbit used in the two previous test 
cases was then numerically integrated for twenty-eight days with this lunar 
gravitational truth model. The satellite orbit was then converted to 
selenographic orbital elements and plotted to determine whether the 
mascons had drastically altered the lunar gravitational field. The plots 
revealed that the mascons affected the satellite when it was in the prime 
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meridian plane, seven and twenty-one days into the orbit propagation. At this 
point, the instantaneous semi-major axis increases noticeably. The lunar 
altitude, however, never drops below 100 km over the integration period. 
The lunar altitude for this orbiter oscillates between 100 - 300 km. When 
numerically integrated without mascons, the satellite orbit exhibits the same 
altitude oscillations, but the semi-major axis does not show the peaks 
observed in the mascon case. The inclusion of mascons does not affect the 
orbital inclination and has a limited effect on the variation of eccentricity 
over the twenty-eight day period. This analysis showed that the addition of 
mascons to the lunar gravitational field did not detrimentally affect the 
satellite orbit. 


Table 6.4-1: "Truth" Model Major Mass Anomalies 


Mascon 

Strength 

xlO ” 6 

Lunar Mass 

Longitude 

Latitude 

# of 
Point 
Masses 

Sea of Rains 

22.8 

328° -350° 

28° - 48° 

u 

Sea of Serenity 

21.5 

8 ° - 25° 

17° - 34° 

12 

Sea of Crises 

9.2 

52° - 63° 

12° - 23° 

jfSji 

Sea of Nectar 

8.4 

27° - 38° 

-17° - -19° 

6 

HUSH 

4.2 

347° -353° 

7° - 13° 

4 

Sea of Moisture 

6.0 

318° -325° 

-29° - -22° 

3 

\msmm 

3.1 

82° - 93° 

ER9 

6 

Subtotal 

75.2 


49 

Difficult to Obs 

40.025 

175° -185° 

-8° - 8° 

19 

Balancing 

27.157 

212° -217° 

-51° - -45° 


Point Mass 

1.114 

217.6° 

-63.4 

B 

| Subtotal 

68.296 


27 
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6.5 Single Orbiter, Earth-Based Doppler Sensing Scheme 

The first sensing scheme evaluated with the lunar gravitational truth 
model was the earth-based Doppler observation of a single lunar orbiter. This 
scheme was selected because it was used to obtain the current lunar gravity 
field models. This scheme's ability or inability to estimate the "true" lunar 
gravitational field may provide an idea of how well current models represent 
the real lunar gravity field. The lunar orbiter was placed in a polar near¬ 
circular orbit which would provide total lunar surface coverage over the 
fourteen day observation period. This orbit attempted to recreate the "gravity 
sensing" satellite's orbit in the Lunar Observer or Scout missions [16, 25, 39], 

Since the acceleration due to gravity is larger for low altitude orbits, its 
disturbances are easier to sense and the lowest possible orbital altitude was 
desired for this mission. Unfortunately this desire conflicts with the desire to 
observe undisturbed motion for at least one lunar period. Very low altitude 
orbits typically require re-boosting to keep them at a safe distance above the 
lunar surface. Re-boost maneuvers, however, disrupt the estimation 
procedure by introducing new forces on the orbiting body which are difficult 
to include in the estimation process. A 100 km altitude polar orbiter was 
numerically integrated for twenty-eight days. Plots of its selenographic orbital 
elements revealed that this orbiter barely remained above the lunar surface 
for the twenty-eight day period. Since this orbiter would require re-boosting 
during the twenty-eight days, the orbital altitude was increased by 100 km. 
The 200 km altitude polar orbiter maintained an altitude of 100-300 km over 
the twenty-eight day period. This satellite orbit was a compromise satisfying 
both the low-altitude and no re-boost requirements. 


6.5.1 Single Orbiter 'Truth" Model Observations 

The single polar 200 km altitude lunar orbiter motion was numerically 
integrated for fourteen days from 16 May 1968 (JD 2,440,001.5) with the 
satellite initial conditions given below. The orbit was propagated using the 
Adams-Moulton numerical integration technique with a step size of 2‘ 13 day. 
The integration step size was relaxed from previous cases to save computer 
disk space, especially when partial derivatives are calculated at each step. 
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Relaxing the integration step size still provided over 725 steps per 2.125504 
hour orbit, more than required for numerical stability and adequate for the 
characterization of an eighth degree and order harmonic fit model. 

Over this fourteen day period, Doppler observations of the lunar 
orbiter were processed every sixty seconds, except during lunar occultations, 
from the DSN Goldstone station using the S-Band frequency (2.115 GHz). 
Transponder frequency translation was not simulated and no horizon 
constraints were imposed. A 14 mHz accuracy was assumed for the Doppler 
observations, roughly equivalent to the quoted 1 mm/sec range rate accuracy 
for sixty second intervals [29]. The simulations of this mission produced 
15,695 Doppler shift observations over the fourteen day period. 

Table 6.5.1-1: Satellite Initial Conditions for 'Truth" Model Observations 


ao 

1938 km (1.295472995 x 10‘ 5 AU) 

eo 

0.05 

io 

90° (103.1048493849350°) 

Do 

90° (304.1996805394721°) 

wo 

90° (72.2394269798987°) 

Mo 

1 ° 


6.5.2 Eighth Degree and Order Fit 

The first estimation run attempted to fit the observations to an 8 X 8 
spherical harmonic expansion. The final iteration of the limb mascon test 
case provided the initial guess of parameters and initial conditions for this 
case's iterations. This saved computer time since both cases were based on 
observations of the same polar orbiter and the limb mascon case had already 
integrated the satellite's motion with partial derivatives. 

Due to the significant difference between the model used to generate 
the "true" observations and the fit model, this case required twenty-five 
iterations to converge upon a solution. The estimated harmonic coefficients 
for the fit model are listed in Table E-3 in Appendix E. The simulation also 
estimated the satellite orbit's initial conditions and Table 6.5.2-1 compares the 
true initial conditions with their estimated values. 
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The observation residuals divided by the assumed measurement 
errors, referred to as the non-dimensional "fit residuals", had a root mean 
square value of 365.954 for this model. This is over thirty times larger than 
the Face Mascon test case and is a further indication of the 8 X 8 spherical 
harmonic expansion's inability to fit the observed satellite's motion. Eighteen 
of the 3,403 parameter correlations were greater in magnitude than 0.90 and 
nine of these were greater than 0.95. Four pairs of correlations (C 55 C 75 , 
S 66 ^ 86 / C 66 C 86 / and S 55 S 75 ) were greater than -0.998 and two pairs (S 54 S 74 and 
C 54 C 74 ) were greater than -0.988. 


Table 65.2-1: 8X8 Single Orbiter Initial Condition Estimates 



"True" 

Fit 

A 

ao 

1938.0 km 

1938.1154 km 

115.4 meters 

eo 

0.05 

0.047811 

-0.002189 

io 

103.1048° 

103.0998° 

-0.0050° 

&0 

304.1997° 

304.1957° 

0.0040° 

(Q-ko)o 

16.4391° 

15.9600° 

-0.4791° 

(Q-ko+M)o 

17.4391° 

17.3912° 

-0.0479° 


6.5.3 Twelfth Degree and Order Fit 

The estimated gravity model's fit residuals show that the 8 x 8 spherical 
harmonic expansion was not a close fit to the observed behavior. Mascons in 
the truth model result in very high frequency gravitational behavior in the 
local area (Section 2.3). Since the spherical harmonic expansion requires 
higher degree and order expansions to model this high frequency behavior, a 
fit to a 12 X 12 spherical harmonic model was attempted. 

For the initial parameter guesses, estimates from the 8 X 8 fit model 
were used and all of the new coefficients (ninth through twelfth degrees) 
were set to zero. The satellite initial condition estimates for the 8 x 8 fit, used 
for the 12 X 12 fit's initial iteration, are listed in Table 6.5.3-1. The same 
"truth" model observations used in the 8 X 8 fit were used for this fit. 
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Since the fit model was increased to twelfth degree and order, the 
numerical integration step size was reduced to 2'^ days. The step size was 
changed to sample the higher degree harmonic effects more often with 1,451 
steps per 2.125504 hour orbit. With the change in step size and dramatically 
increased number of estimation parameters and parameter partials, the 
numerical integration required approximately 23 hours to propagate the 
satellite's motion and partial derivatives for fourteen days. 

Table 6.5.3-1: 12 X 12 Fit Initial Guesses for Satellite Initial Conditions 


ao 

1.29555011876471 x 10' 5 AU 

eo 

0.0478107828679364 

io 

103.099831198443° 

Do 

304.195650955435° 

(Q. + co)o 

15.9599887633931° 

(Q + co + M)o 

17.3911644150588° 


On the fitting process' second iteration several of the parameter 
adjustments were an order of magnitude larger than their previous estimate. 
These estimates were used to propagate the next iteration's orbit, resulting in 
increased residuals and even larger parameter adjustments. By the fifth 
iteration, the normal equations could not be solved because 152 of the 171 
diagonal elements in the coefficient matrix (A) were negative. Since the 
normal equations could not be inverted, the fitting process was abandoned. 


Table 6.5.3-2: 12 X 12 Spherical Harmonic Fit Progression 


Iteration 

Pre-adjustment 
RMS Residual 

Predicted RMS 

Residual 

1 

365.954 

112.113 

2 

10,433.0 

99.5596 

3 

42,590.9 

192.212 

4 

106,490 

1,999.54 

5 

760,314,000 

16,484,700 


Table 6.5.3-2 shows how the parameter adjustments kept leading the 
process further and further from a solution. The fitting process started with 
observation residuals divided by the assumed measurement errors with a 
root mean square value of 366. After solving the normal equations, PEP 
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predicted a non-dimensional rms residual of 112. The predicted residual is 
the computed pre-adjustment residual less the sum of the observable partial 
derivatives times the parameter adjustments. Instead of reducing the 
residuals by a third, the new parameter estimates produced residuals with a 
non-dimensional rms almost thirty times greater than in the previous 
iteration! Unfortunately instead of making it easier to fit to the observed 
lunar orbiter motion, increasing the degree and order of the spherical 
harmonic model made it more difficult to solve the normal equations and 
converge upon a fit solution. This may have occurred because of the limited 
observability of lunar far-side motion and high parameter correlations. 


6.6 Tenth Degree Harmonic Truth and Fit Model Test Case 

After the previous convergence difficulties, a test was run to estimate a 
10 X 10 spherical harmonic expansion for observations generated with the 
same degree and order truth model. Since the previous case was the first 
attempt to fit to harmonic expansions higher than eighth degree and order, 
this test would verify whether or not the difficulties encountered were due to 
the estimation program or the spherical harmonic expansion's ability to fit to 
the observed behavior. 

The same lunar polar 200 km altitude near-circular orbit used in 
previous runs was used to generate Doppler "truth" observations with 1 
mm/sec accuracy. This orbit, again integrated for fourteen days from 16 May 
1968, used the orbital initial conditions listed in Table 6.2-1 and a 10 X 10 
spherical harmonic lunar gravity truth model. This truth model used the 
1980 Bills and Ferrari coefficients [12] for the initial five degrees (Table 6-2.2). 
Gravitational coefficients developed at JPL by Alex Konopliv for his 50 X 50 
spherical harmonic expansion were used for the sixth through tenth degree 
harmonic coefficients [27], because the Bills and Ferrari coefficients [12] were 
not manually entered into the PEP input stream. The coefficients used are 
listed in Table E-4 in Appendix E. 

After generating the "truth" model observations, the satellite initial 
conditions and one hundred and seventeen lunar harmonic coefficients were 
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perturbed to provide the first guess for the estimation process. As in the 
previous fit to a pure harmonic truth test case, each gravitational harmonic 
coefficient's absolute value was increased by 1.0 x 10' 7 and the orbital initial 
conditions were perturbed by the values given in Table 6.2-3. The satellite 
orbit was then numerically integrated with partial derivatives. 

In the first iteration the pre-adjustment observations residuals divided 
by the assumed measurement error had a root mean square value of 715,121! 
As a result, the calculated parameter adjustments were orders of magnitude 
larger than the guessed parameters. After the first iteration, the lunar orbiter 
eccentricity grew from 0.05001 to 0.24325 and the semi-major axis was 
increased by 532 km. These initial conditions were then propagated to create 
the next iteration's observations. In the second iteration, the normal 
equations could not be solved because 51 of the 123 diagonal elements in the 
coefficient matrix (A) were negative. Once again since the normal equations 
could not be inverted, the estimation process was abandoned. Rather than 
converging upon a solution. Table 6.6-1 demonstrates how quickly the process 
diverged. 


Table 6.6-1: 10 X 10 Spherical Harmonic Fit Progression 


Iteration 

Pre-adjustment 
RMS Residual 

Predicted RMS 

Residual 

1 

715,121 

248,467 

2 

3.20623 x 10 10 

2.11726 x 10 7 


Since the previous test case (Section 6.2) encountered difficulty when 
an incomplete set of harmonic coefficients were used to generate partial 
derivatives of the satellite motion, this cases' observation partial derivatives 
were checked using a finite difference method (Section 3.4.4). This check 
verified that the observation partials were correct. 

Since the estimation software was operating properly, the gravitational 
model's sensitivity to initial guesses seemed to cause the convergence 
difficulties. Apparently the initial guesses resulted in theoretical observation 
values so different from the observed behavior that correct parameter 
adjustments could not be determined. Since the estimation of higher degree 
and order expansion fit models is extremely sensitive to the initial parameter 
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guesses, fit models above eighth degree and order were not attempted for the 
remainder of this thesis' simulations. 

Smaller parameter adjustments from iteration to iteration could help 
the process converge. Incorporating a priori estimates of the lunar 
gravitational field and their uncertainties would reduce the size of the 
parameter adjustments (see Section 5.1.4). Smaller parameter adjustments 
could also be used by under-weighting the calculated adjustments, i.e. only 
adjusting the initial guesses by 2/3 of the calculated adjustment on the first 
iteration, etc. Convergence could also be aided by first estimating an 8 X 8, 
then a 9 X 9, and then a 10 X 10 degree and order fit model to the observed 
behavior, using the estimates from one model as initial guesses for the next 
higher degree model. 


6.7 Dual Orbiter, Bent Pipe Doppler Sensing Scheme 

The next sensing scheme evaluated with the lunar gravitational 
"truth" model featured two lunar satellites and the simulation of a bent pipe 
Doppler sensing scheme between earth-based sites and the two satellites. This 
sensing scheme was selected since it is one of the schemes being considered by 
NASA for future lunar gravitational sensing missions. 

This sensing scheme uses a low altitude circular polar "gravity 
sensing" satellite. A coplanar elliptical "viewing" satellite makes lunar far- 
side observations of the "gravity sensing" satellite's motion [16, 39, 40]. The 
polar 200 km altitude near-circular lunar orbiter used in previous 
simulations is once again the "gravity sensing" satellite for this dual orbiter 
sensing scheme. The satellite initial conditions and numerical integration 
parameters for this simulation were the same as for the single orbiter, earth- 
based Doppler observation sensing scheme and are given in Section 6.5. The 
"viewing" satellite was placed in a 450 km X 7,000 km altitude elliptical orbit 
with a 10.06 hour period. This satellite was given the same initial inclination 
(i) and longitude of the ascending node (£2) as the circular satellite, placing 
them in the same orbital plane, orbiting the moon in the same direction. 
With this selection of ascending node, the far side of the moon rotates 
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underneath the elliptical orbit's apolune during the first fourteen days of orbit 
propagation. This maximizes the satellite's ability to view the low altitude 
orbiter's far-side motion. 

The argument of perilune was selected to skew the ellipse and bring 
the orbit's apolune position out of the earth occultation zone (Figure 1.4-1). 
Keeping the satellite's apolune passage out of the occultation zone increases 
the time earth-based sites can view the elliptical satellite for bent-pipe 
measurements as well as data transfers. The long term behavior of the 
osculating orbital elements was studied to determine this skew direction. 
Formulas for the doubly averaged effect of the earth upon a lunar orbit 
revealed that if the initial argument of perilune, to, is in the second or fourth 
quadrants, the polar orbit's eccentricity will decrease and co will drift to the 
first or third quadrant [9]. Based upon this analysis and the desire to keep the 
apolune over the backside of the moon, a perilune angle in the fourth 
quadrant was selected. The initial selenographic orbital elements for this 
elliptical satellite are given in Table 6.7-1, with the angle values in 
parentheses referred to the mean equinox and equator of the earth of 1950.0. 
These initial conditions were numerically integrated for fourteen days from 
16 May 1968 (JD 2,440,001.5) using the Nordsieck variable step size integration 
technique (Section 3.1). Since the elliptical orbit was newly propagated, its 
numerical integration file was converted to selenographic orbital elements. 
Plots of these elements revealed that the orbit was stable and, as predicted, 
the orbital eccentricity did decrease over the fourteen day period. 

Table 6.7-1: Elliptical Satellite Initial Conditions for 'Truth" Model Observations 


ao 

5463 km (3.651789462 x 10' 5 AU) 

eo 

0.05994874611 

io 

90° (103.1048493849350°) 

fio 

90° (304.1996805394721°) 

0)0 

315° (297.2394269798987°) 

Mo 

1 ° 
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6.7.1 Dual Orbiter 'Truth" Model Observations 

The coherent bent pipe observation method (Section 4.2.3) was 
simulated over the fourteen days the satellites' orbits were propagated. One 
two-way coherent Doppler loop processed observations between the DSN 
Goldstone site and the elliptical satellite. Simultaneously, the elliptical 
satellite generated two-way coherent Doppler observations of the low-altitude 
orbiter. Although the coherent bent pipe scheme would be interrupted if any 
of the links were occulted, only the individual links in this simulation were 
interrupted by lunar occultations. In addition to the bent pipe observation 
simulation, the DSN Goldstone site generated Doppler observations of the 
near-circular polar satellite during its near-side passes. All of these 
observations were simulated using the 2.115 GHz S-band frequency with 
approximately 1 mm/sec range rate accuracy (14 mHz) [29]. 

Simulating this sensing scheme produced 15,695 earth-based Doppler 
observations of the circular satellite. Because the elliptical orbit was skewed, 
19,699 Doppler observations were generated between it and the DSN station. 
Over this same period, 10,680 Doppler observations were simulated between 
the two satellites. Auxiliary software analyzed these observation series and 
their occultation periods to evaluate the far-side lunar coverage. Over the 
fourteen days, the low altitude polar orbiter passed behind the moon 114 
times. For most of these occultations, the elliptical satellite viewed this far- 
side passage. The moon blocked the line-of-sight between the two satellites 
on 49 of these 114 occasions and this blockage usually only affected a portion 
of the passage. There were only 15 cases in which the line-of-sight to the 
elliptical satellite was blocked for the entire far-side passage. Since these gaps 
in observation coverage did not occur for sequential far-side passes, the 46,074 
observations should provide excellent visibility into the circular orbiter's 
motion on the lunar far side. 


6.7.2 Eighth Degree and Order Fit 

These observations, generated with the lunar gravitational "truth" 
model, were then fit to an 8 X 8 spherical harmonic expansion representing 
the moon's gravitational field. The initial conditions for the two satellites 
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were perturbed from their true values for the initial guess in the fitting 
process. Table 6.2-3 lists the perturbations which were applied to both the 
circular and elliptical orbits' initial conditions. The true values of the 5 X 5 
spherical harmonic expansion used in the truth model were used as initial 
guesses for the first five degrees of gravitational coefficients (Table 6.2-2). 
Gravitational coefficients from Alex Konopliv's 50 X 50 spherical harmonic 
model of the moon [27] were used for the sixth, seventh, and eighth degree 
coefficient initial guesses. Table E-4 in Appendix E lists these coefficients. 

The iteration fit did not converge upon a solution in its normal sense. 
With the increased observability of lunar far-side motion, the parameter 
estimation routine drove the non-dimensional residuals down to a root 
mean square of 386.462 from an initial 23,763.6 in the first four iterations. 
Because of the difference between the "truth" and fit models, the estimation 
routine could not reduce the non-dimensional residual rms below 360 as 
shown in Table 6.7.2-1. Once again, the predicted residual is the computed 
pre-adjustment residual minus the sum of the observable partial derivatives 
times the parameter adjustments. 

Table 6.7.2-1: 8x8 Spherical Harmonic Fit Progression, 


Iteration 

Pre-adjustment 
RMS Residual 

Predicted RMS 

Residual 

1 

23,763.6 

817.104 

2 

2,983.65 

437.085 

3 

1,746.75 

383.599 

4 

386.462 

371.811 

5 

375.229 

366.587 

6 

368.213 

363.927 

; 

- . - - 

• 

15 

361.210 

361.015 

; 

; 

; 

38 

361.099 

360.939 


Although the estimation routine had minimized the fit residuals, its 
convergence criteria is based on the ratio of the parameter adjustments to 
their uncertainties with the assumed measurement errors, and these values 
remained above the convergence limit of 0.01. Since the convergence 
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criterion did not signal convergence, the estimation iterations continued ad 
infinitum, slowly reducing the residual rms by hundredths and thousandths. 
The process was stopped after thirty-eight iterations. 

The estimated harmonic coefficients for the fit model are listed in 
Table E-5 in Appendix E. The simulation also estimated the orbital initial 
conditions for both of the satellites in the sensing scheme. Tables 6.7.2-2 and 
6.7.2-3 compare the true initial conditions with their estimated values. 

The estimated model's fit residuals had a root mean square value of 
361.099 as listed in Table 6.7.2-1. This is smaller than the root mean square of 
the residuals achieved with the fit to a single lunar polar orbiter (365.954). 
The high fit residuals are an indication of the 8X8 spherical harmonic 
expansion's inability to fit the observed motion due to the "true" 
gravitational field. This dual satellite observation method has significantly 
reduced some of the high parameter correlations, although the most 
significant ones still remain. For this fit case fifteen of the 3,916 parameter 
correlations were greater in magnitude than 0.90 and ten of these were greater 
than 0.95. As with the previous 8x8 spherical harmonic fit, the four highest 
correlations were among the C 55 C 75 , S 55 S 75 / S 6 6S 86 , and C 6 6C 86 pairs. 
Although the correlations between the fifth degree, fifth order and seventh 
degree, fifth order terms were still greater than -0.998, the other correlations 
were now reduced to -0.995. The S 54 S 74 and C 54 C 74 pairs were again the next 
highest correlated (-0.988). 


Table 6.7.2-2: Circular Orbiter Initial Condition Estimates 



"True" 

Fit 

A 

ao 

1938.0 km 

1937.9401 km 

-59.9 meters 

eo 

0.05 

0.049677 

-0.000323 

io 

103.1048° 

103.1033° 

-0.0015° 

Qo 

304.1997° 

304.1994° 

-0.0003° 

(ft-KO)o 

16.4391° 

16.2649° 

-0.1742° 

(Q-KCH-M)o 

17.4391° 

17.4469° 

0.0078° 
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Table 6.7.2-3: Elliptical Orbiter Initial Condition Estimates 



"True" 

Fit 

A 

ao 

5463.0 km 

5463.6291 km 

629.1 m 

eo 

0.599487 

0.599573 

0.000086 

k 

103.1048° 

103.1068° 

0 .0020° 

Qo 

304.1997° 

304.1992° 

0.0005° 

0)0 

297.2394° 

297.2542° 

0.0148° 

Mo 

1 .0° 

0.9997° 

-0.0003° 
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Chapter Seven 


Lunar Gravity Field 
Estimation Analyses 


7.1 Analysis Techniques 

After converging upon a least squares maximum likelihood estimate 
of gravitational coefficients, the estimated lunar gravity fields were analyzed 
to compare their ability to estimate the "true" gravity field. Since spaceborne 
navigation depends on accurately modeling the forces acting on a spacecraft, 
these analyses focus on the effects of modeling errors on lunar navigation. 

The root mean square of the observation residuals does not provide an 
adequate analysis of the global lunar gravity errors due to mismodeling, since 
it only considers the areas of the gravity field where observations were made. 
This underweights the far side of the moon where fewer measurements are 
taken. In Section 7.2 the global radial accelerations for the estimated fit lunar 
gravity fields are compared to the "true" lunar radial accelerations. The limb 
mascon and face mascon estimated fit models are analyzed and the two 8x8 
spherical harmonic expansion fit models are compared to the lunar 
gravitational "truth" model developed in Section 6.4. 

Next, two different lunar spacecraft mission phases are simulated to 
evaluate the two estimated fits to the lunar gravitational "truth" model. 
These analyses show how state errors grow using the estimated gravity field 
model and are intended to simulate the real-world consequences of planning 
and executing lunar missions with a mismodeled gravity field. 

In Section 7.3 the state uncertainties for a low inclination, low altitude 
spacecraft orbit are predicted one orbit ahead using the covariance analysis 
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described in Section 5.2. Stations tracking a lunar satellite would use this 
analytical method during far-side passage to predict the spacecraft's state 
uncertainties upon reemergence from the backside. Since this state 
uncertainty prediction is based upon a mismodeled lunar gravity field, the 
predicted uncertainties will not be exact. The "true" state is calculated by 
numerically integrating the spacecraft orbit using the lunar gravity "truth" 
model. The difference between the spacecraft's predicted and "true" states is 
the error. Comparing the "true" state with the predicted states reveals the 
accuracy of the estimated lunar gravity field. 

A lunar landing deorbit maneuver is simulated in Section 7.4. In this 
analysis the deorbit trajectory from a low inclination circular parking orbit 
was determined from the estimated lunar gravity field models. The mission 
was then propagated using the "true" lunar gravity field. Lunar gravity 
mismodeling resulted in a spacecraft position error by the time the spacecraft 
reached the lunar altitude for Powered Descent Initiation (PDI), typically 
about 18 kilometers. In a real lunar mission, the spacecraft would be forced to 
burn extra propellant in a suboptimal descent trajectory to recover from these 
errors. 


Finally, the high gravitational coefficient correlations encountered in 
the estimation process are analyzed in Section 7.5. Three different 
measurement types and orbital orientations are simulated to determine their 
effect on the correlations. Reducing the coefficient correlations by 
introducing new observations may permit the process to converge more 
quickly and estimate a more accurate gravitational model. Additionally, the 
highest correlations are studied in an attempt to find ways to allow them to be 
estimated independently. 


7.2 Global Lunar Radial Acceleration Analysis 

Although the estimation process provides residual statistics, they only 
provide an idea of how well the estimate fits the observations since the 
observations are not available over the entire lunar surface. To analyze the 
fit globally a software program was written to calculate the radial accelerations 
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for a gravity field model everywhere on a sphere of constant radius. Only the 
radial accelerations due to gravity harmonics and/or mascons are calculated 
in this program. These radial accelerations were calculated for fit and truth 
models for all of the cases estimated in Chapter Six at a lunar altitude of 100 
km for grid points spaced every four degrees of selenographic latitude and 
longitude. Unfortunately, since the commercial graphics package was 
erroneously plotting lines of constant radial acceleration, contour plots 
comparing the fit and “truth" models are not currently available. 

Based on the calculated radial accelerations, the fit and “truth" models 
were compared statistically. The radial acceleration errors between the 
"truth" and fit models were calculated for each grid point. The root mean 
square error between the two models was then calculated. Because of the 
even spacing of grid points in selenographic latitude and longitude (chosen 
for rectangular contour plots), this analysis weights radial acceleration errors 
more heavily in the polar regions. For the mascon test case this bias does not 
favor either scenario since the disturbing mascons were placed on the 
equator. For the two 8x8 estimated fit models, the dual orbiter bent pipe 
estimated fit might have an advantage since the "viewing" satellite observed 
the “sensing" satellite as it passed over the lunar polar regions. Both models, 
however, were trying to estimate the same lunar gravity truth model, so a 
comparison of the rss radial acceleration errors is still valid. 


7.2.1 Limb / Face Mascon Analysis 

Based on the global lunar radial acceleration errors, the face mascon 
case produced a slightly better estimated gravity field. This case had an rms 
radial acceleration error at 100 km altitude of 64.7804 milligals. The estimated 
fit to the limb mascon case had radial acceleration errors with an rms of 
66.6636 milligals. This global analysis of the radial accelerations errors reveals 
that a slightly better lunar gravity field was estimated when the disturbing 
mascons were aligned with the sensing line-of-sight. 

This test case attempted to determine whether the mascon's location 
and the sensing geometry affected the estimated gravity field. The non- 
dimensional fit residuals' rms contradicted the global radial acceleration 
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analysis. The limb mascon fit model (7.37) resulted in a better fit than the face 
mascon fit model (11.80). These statistics could be the result of the near-side 
face mascon causing larger line-of-sight accelerations for the orbiting satellite 
than the limb mascons produced. These statistics could also result from the 
measurement's ability to sense the near-side mascon's disturbance or inability 
to recognize any disturbance from the limb mascons. Unfortunately, this 
analysis does not indicate which estimated field is better at modeling the 
radial accelerations in the local vicinity of the mascons. Contour plots of 
these radial accelerations will demonstrate the fit models' ability to locally 
model these mascon disturbances. 


7.2.2 Analysis of the Eighth Degree and Order Fits 

Based on the global lunar radial accelerations errors, the dual orbiter 
bent pipe sensing scheme produced the best estimated gravity field. This case 
had an rms radial acceleration error at 100 km altitude of 263.015 milligals. 
The estimated fit to the single orbiter, earth-based sensing scheme had radial 
acceleration errors with an rms of 320.073 milligals. The disparity between 
the two rms values provides a strong indication of the advantage of including 
lunar far-side observations in the estimation process. These errors also give 
an indication of how much more difficult it was to estimate a spherical 
harmonic expansion in the presence of 79 point mass disturbances than it was 
in the mascon test case with just two point masses. 


7.3 Single Orbit State Uncertainty Prediction Analysis 

For the next analysis, a low inclination, low altitude lunar satellite 
orbit was used to analyze the position and velocity errors between the "true" 
and estimated gravity fields. The lunar gravity field truth model developed 
in Section 6.4 and the two 8X8 spherical harmonic expansion fit models 
derived from the two different sensing schemes were used for this analysis. 
A 15° inclination, 100 km lunar altitude satellite orbit was selected for this 
analysis. The estimated 8X8 spherical harmonic fit models and their 
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coefficient uncertainties were first used to predict the state uncertainties for 
this orbit. This covariance analysis was performed with both the gravity 
harmonic coefficient covariance matrix (A' 1 ) and this matrix multiplied by 
the fit residuals' rms (X). Fit residuals, once again, are the observation 
residuals divided by the assumed measurement errors. This evaluation orbit 
was then numerically integrated using the "true" lunar gravity field. The 
state errors between the "true" orbits and those predicted from the two 8 X 8 
estimated spherical harmonic models provided a further measure of the two 
sensing schemes' capabilities and limitations. Furthermore, the predicted 
uncertainties were compared to the errors to understand the uncertainty 
prediction accuracy with gravity field modeling errors. 

The equations of motion and equations of the partial derivatives of 
motion with respect to the gravity harmonic coefficients were numerically 
integrated for one orbit (117.85 minutes) for the two estimated fit models 
from 16 May 1968 (JD 2,440,001.5 0 h CT) with the satellite selenographic initial 
osculating orbital elements listed in Table 7.3-1 (1950.0 angles in parentheses). 
Based on these initial conditions, the lunar orbiter has emerged from behind 
the far side of the moon and is in the middle of a near-side pass. In the 
numerical integration of the "truth" model, partial derivatives were not 
calculated. Auxiliary software programs then used the numerical integration 
output files for the three cases as well as the covariance files for the two 
estimated fit models to perform the error analysis. 

Table 7.3-1: Low Inclination, Low Altitude Evaluation Orbit 
Satellite Initial Conditions 


ao 

1838 km (1.22862712 x 10* 5 AU) 

eo 

0.01 

io 

15° (30.7285027175189°) 

Q) 

210° (26.9029402084826°) 

too 

180° (222.6750569007985°) 

Mo 

o 

in 

H 

CO 
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7.3.1 State Uncertainty Prediction 

The state uncertainty prediction based on PEP's gravity harmonic 
coefficient covariance matrix (A -1 ) derived from the assumed measurement 
errors was overly optimistic for both cases since the fit residuals for the two 
estimated runs were so high (365.954 and 361.099 rms). This method 
predicted root sum squared (rss) state uncertainties less than 500 meters in 
position and one half of a meter per second in velocity over the orbit for the 
single orbiter estimated fit case. For the dual orbiter estimated fit model this 
method predicted rss uncertainties of less than 250 meters in position and one 
quarter of a meter per second in velocity over the same orbit. These 
predictions are extremely optimistic, especially when compared to the orbital 
errors presented in Section 7.3.2. 

The prediction based on the estimated fit's covariance matrix £ from 
Equation (5.2-2) provided a more realistic assessment of the state 
uncertainties. This method predicted rss state uncertainties as high as 4 
kilometers in position and 4 meters per second in velocity for the single 
orbiter estimated fit case. These predicted state uncertainties, transformed to 
the local vertical, local horizontal coordinates, are plotted over time for this 
lunar orbit in Figure 7.3.1-1. For the dual body estimated fit model, this 
analysis predicted rss uncertainties of approximately 750 meters in position 
and 0.6 meters per second in velocity. These local vertical, local horizontal 
predicted state uncertainties are plotted over time in Figure 7.3.1-2. These 
state uncertainty plots show that the largest position uncertainties are in the 
range direction and the largest velocity uncertainties are in the vertical 
direction as might be expected for the radial acceleration perturbations caused 
by the mascons. 

This state uncertainty prediction analysis demonstrates the importance 
of including lunar far-side observations in the estimation process. The 
uncertainties predicted for the gravity model estimated from the dual orbiter 
sensing scheme were approximately one fourth of those estimated from the 
single orbiter sensing scheme. The dual orbiter cases' uncertainties were also 
very low while the lunar satellite passed across the lunar near side (first forty 
minutes of numerical integration). For both cases, the uncertainties grow as 
the orbiter passes behind the far side of the moon, reaching local maximums 
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for velocity in the vertical direction and position in the down range direction 
in the middle of the far-side pass. These uncertainties subside when the 
spacecraft reemerges from the far side of the moon eighty-eight minutes into 
the orbit. From this point the uncertainties increase during the final lunar 
near-side pass. 


7.3.2 State Errors Between "True" and Estimated Gravity Models 

The errors between the "true" and estimated fit lunar gravity models 
were determined from a direct comparison of the two numerical integration 
files. These state errors were then transformed to local vertical, local 
horizontal coordinates. The local vertical, local horizontal state errors 
between the orbits generated for the single orbiter estimated lunar gravity 
field and the "true" gravity field are plotted in Figure 7.3.2-1. The same errors 
resulting from the difference between the dual orbiter estimated gravity field 
and the "true" gravity field are plotted versus time in Figure 7. 3 . 2 - 2 . 

The errors for the single orbiter estimated gravity field stay relatively 
flat as the orbiter crosses the near side of the moon. Forty minutes into the 
numerical integration the orbiter passes behind the moon. Midway through 
the far-side pass there is a rapid growth of velocity errors. Shortly thereafter 
these velocity errors manifest themselves as position errors as large as 15 
kilometers in the range and 10 kilometers in the cross track directions. These 
errors demonstrate this estimated fit model's inability to estimate the far-side 
gravity field. 

The errors for the dual orbiter estimated gravity field are much lower 
than in the previous case. When plotted on the same scales as the errors in 
the first case, the dual orbiter case's velocity errors stay in a narrow band 
around zero meters per second. The dual orbiter position errors are also 
much more reasonable for the cross track and vertical directions. The dual 
orbiter estimated gravity model still results in large range errors over the 
single lunar orbit. Since the dual orbiter estimated lunar gravity model is 
based on observations over the entire lunar surface, the position and velocity 
errors do not manifest themselves at a specific point in the lunar orbit as was 
the previous case. In this case the cross track and vertical position errors and 
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Single Orbit Spacecraft State Uncertainty 


Predicted State Uncertainty for Single Orbiter Sensing Fit 
Inclination = 15 deg, Eccentricity = 0.01, Altitude = 100km 




Figure 7.3.1-1: Single Orbiter Estimated Gravity Field State Uncertainties 
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Single Orbit Spacecraft State Uncertainty 

Predicted State Uncertainty for Dual Orbiter Sensing Fit 
Inclination = 15 deg, Eccentricity = 0.01, Altitude = 100km 




Figure 7.3.1-2: Dual Orbiter Estimated Gravity Field State Uncertainties 
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Single Orbit Spacecraft State Errors 

Errors Between Single Orbiter Sensing Scheme Fit and Truth Model 
Inclination = 15 deg, Eccentricity = 0.01, Altitude = 100km 




Figure 7.3.2-1: Single Orbiter Estimated Gravity Field State Errors 
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Single Orbit Spacecraft State Errors 

Errors Between Dual Orbiter Sensing Scheme Fit and Truth Model 
Inclination = 15 deg, Eccentricity = 0.01, Altitude = 100km 




Figure 73.2-2: Dual Orbiter Estimated Gravity Field State Errors 
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the range and cross track velocity errors oscillate about zero for the single 
lunar orbit. The vertical velocity errors and resulting down range position 
errors, however, grow without oscillating about zero for the single orbit 
propagation. 

For both cases the state uncertainty prediction was very optimistic, 
even after including the rms of the fit residuals. This discrepancy results 
from the difference between the structure of the two lunar gravity models 
("truth" and fit). The rss of the predicted position uncertainty was an order of 
magnitude less than the "true" position errors. The predicted uncertainty in 
velocity was only slightly better. Additionally, the covariance analysis does 
not successfully predict cross track position uncertainties. This is especially 
apparent for the single orbiter estimated gravity field case. Figure 7.3.2-1 
shows significant cross track state errors and no notable cross track 
uncertainties are predicted in Figure 7.3.1-1. 

When analyzing the proposed Lunar Observer mission with different 
truth and fit models, Alex Konopliv also noted that a covariance uncertainty 
analysis was unbelievable because it was "overly optimistic for all cases" [25]. 
This deficiency in the covariance analysis is an indication that the system 
dynamic model needs to include some process noise to account for the gravity 
field mismodeling. A Kalman filtering or maximum likelihood system 
identification technique could include this process noise as it propagated the 
satellite equations of motion. 


7.4 Lunar Deorbit Maneuver Error Analysis 

Finally, a satellite lunar landing from a low inclination circular orbit 
was simulated for both estimated gravity fields. The deorbit burn for this 
lunar maneuver was determined based on numerical orbit integrations using 
the estimated lunar gravity fields. This maneuver was then executed in the 
"true" lunar gravity field, and the position errors at PDI were used to evaluate 
the two estimated gravity field models' ability to plan future lunar missions. 
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This analysis assumed that the spacecraft was resupplying a lunar base 
on the near side of the moon and had been inserted into a circular, low 
inclination orbit. The transfer from this orbit to the lunar surface was 
planned using a Hohmann transfer. The circular parking orbit was 
numerically integrated from the near side of the moon to the middle of the 
far side of the moon. At this point, a deorbit burn placed the spacecraft in an 
elliptical transfer orbit. The spacecraft's pre-planned Powered Descent 
Initiation (PDI) location coincided with the elliptical transfer orbits' perilune. 
From this point, the spacecraft would begin its powered descent to the lunar 
surface and rendezvous with the lunar base. The powered descent portion of 
the landing mission was not simulated in this analysis. 



Figure 7.4-1: Lunar Deorbit Mission 


Earth 


The mission simulation began by integrating the spacecraft's circular 
orbit starting on 16 May 1968 (JD 2,440,001.5 Oh) with the satellite selenographic 
initial conditions listed in Table 7.4-1 (1950.0 angles in parentheses). Based on 
the satellite's state at M=180° in the parking orbit, an initial guess for the Av 
was calculated for a Keplerian transfer orbit with the PDI perilune altitude, 
the lunar surface for this analysis. The elliptical transfer orbit was then 
numerically integrated and the guessed Av updated based upon the PDI 
position error at perilune. These iterations were continued until acceptable 
target PDI positions were obtained. Since this maneuver was planned for the 
two different estimated lunar gravity fields, two different Av's were calculated 
and subsequently executed. 
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For the single body estimated gravity field case, the calculated deorbit 
burn occurred over 0.0611° south latitude and 179.9215° west longitude 
approximately 63.985 minutes into the numerical integration. The transfer 
Av had a magnitude of 2.535 x 10' 5 AU/day (PEP's units) with a projected PDI 
at 0.040° south latitude and 1.951° west longitude. This burn was then 
executed by the spacecraft in the "true" lunar gravity field. At 63.985 minutes 
into the numerical orbit integration the previously calculated Av was 
subtracted from the satellite's inertial velocity to simulate the deorbit burn. 
The numerical integration then proceeded from this new state. This 
maneuver resulted in a PDI at 0.766° south latitude and 9.903° west longitude, 
a 7.952° error in longitude and 0.725° error in latitude! 

Table 7.4-1: Low Inclination, Lunar Landing Parking Orbit 
Satellite Initial Conditions 


ao 

1938 km (1.295472995 x 10' 5 AU) 

eo 

0.0005 

io 

5° 

(18.2129222440553°) 

Do 

0° 

(349.1664497093493°) 

too 

0° 

(226.5049507382479°) 

M 0 

0° 


The deorbit was then planned with the dual orbiter estimated lunar 
gravity field. The deorbit burn was again scheduled for 63.985 minutes into 
the orbit (from the satellite initial conditions in Table 7.4-1) at a selenographic 
latitude of -0.0358° and longitude of -179.9851°. The planned magnitude of 
the Av was 2.5225 x 10* 5 AU/day with a projected PDI of 0.4339° south latitude 
and 5.8691° west longitude. This maneuver was then executed with the 
"true" lunar gravity field. This simulation resulted in a PDI at 0.2503° south 
latitude and 4.008° west longitude. The PDI errors in this case were now 
1.8613° in longitude and 0.1835° in latitude, a significant reduction in the PDI 
position errors and further evidence of the importance of including lunar far- 
side observations in lunar gravity field estimations. 

On the lunar surface these errors in latitude and longitude would 
result in a position error of 242 kilometers in the single orbiter estimated 8X8 
field case and 56 kilometers for the dual orbiter estimated 8x8 field. If the 
landing spacecraft began powered descent this far off target it should still be 
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able to land close to the lunar base, as long as it is close enough to receive a 
radio beacon signal. The suboptimal descent trajectory required to 
compensate for these errors, especially in the first case, would seriously cut 
into any propellant margins, perhaps jeopardizing future mission plans. 

The magnitude of the navigation errors might be reduced if higher 
than 8x8 spherical harmonic or other gravity model fits were estimated. For 
higher degree and order fits, it would be interesting to note whether the 
above factor of 4.3 improvement in navigation accuracy still held when far- 
side satellite-to-satellite measurements were added to near-side earth-based 
measurements. 


7.5 Gravity Coefficient Parameter Correlation Analysis 

Each of the estimation runs in Chapter Six produced extremely high 
correlations for some of the gravitational coefficients. These high 
correlations inhibited the estimation routine's ability to converge quickly and 
provide an accurate estimate of the gravitational parameters. Different 
measurement types and orbital orientations were analyzed to see if any of 
these parameter correlations could be broken. For this analysis, truth 
model observations were generated for the new observation methods. Based 
on each new set of observations, the lunar gravity field was estimated for a 
single iteration of the process outlined in Figure 6.1-1. The gravitational 
parameter correlations were then obtained from the single iteration s 
covariance matrix. 

For the first correlation analysis, the elliptical "viewing" orbit in the 
dual orbiter, bent pipe observation scheme was placed in an orthogonal, 
rather than coplanar, orbit plane. The "sensing" orbit was left in its polar 
orbit to provide full lunar surface coverage over the fourteen day observation 
period, so the elliptical satellite was placed in an equatorial orbit. From this 
orbital geometry, bent pipe observations identical to the coplanar case were 
simulated for the estimation process. From this iteration, sixteen of the 3,926 
correlations were greater in magnitude than 0.90. Fourteen of these were 
greater than 0.95. The C 55 C 75 , S 55 S 75 pairs had the highest correlation of 
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-0.999. The 054074 , S54S74 correlation pairs were the second highest at -0.996. 
The C 66086 and §66 §86 pairs also had correlations greater than -0.993 

The second method investigated an interferometric observation 
method to determine its impact on the high gravitational parameter 
correlations. For this scenario, NASA DSN stations made milli-arc second 
long baseline interferometer measurements of the near-circular 200 km 
altitude lunar orbiter in addition to the 3 meter range and 1 mm/sec Doppler 
observations. Because of difficulties encountered simulating PEP's internal 
interferometry observations, these interferometer measurements were 
simulated by azimuth and elevation angle observations for a single DSN site 
(Goldstone) with the interferometric accuracy. Since the interferometer 
observations only provide a planar angular measurement, the fitting process 
was executed ignoring the elevation angle observations. The fitting process 
was then repeated using both of the angle measurements to see what impact 
precise earth-based three-dimensional angular observations combined with 
range and Doppler observations would have on the gravitational parameter 
coefficients. 

For both of these fit cases (azimuth only/azimuth and elevation), 
twenty two of the 3,403 correlations were greater in magnitude than 0.90 and 
fifteen of these were greater than 0.95. The parameter correlations were 
almost identical for the two cases, with the highest parameter correlations 
differing only in the fourth or fifth decimal place. The S55S75 and C55C75 
pairs were the most highly correlated (-0.999). The C 65 Css, S 6 5S 8 5, C54C74, 
S54S74 parameter correlations were the second highest group. For the 
interferometric case the 66 , 86 correlations are still very high (> -0.990). 

For the final case, the Goddard Space Flight Center's (GSFC) proposed 
co-orbital laser ranging and Doppler lunar gravity field sensing scheme was 
simulated (Figure 1.4-2). The main satellite was placed in the polar near- 
circular 200 km altitude orbit which has proven so useful throughout this 
thesis. The co-orbiting subsatellite was placed in the identical orbit but with 
an initial mean anomaly nine degrees ahead of the main satellite. Laser 
ranging and sixty second Doppler observations were processed for the entire 
fourteen day orbit with 1 mm laser ranging and 1 mm/sec Doppler accuracies. 
These accuracies were based on GSFC briefed capabilities [ 2 ]. Additionally this 
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method used two-way coherent Doppler observations between the DSN 
Goldstone site and the main satellite using the S-Band (2.115 GHz) frequency 
with 1 mm/sec accuracy for sixty second count intervals. 

For this case, the parameter correlations were different from all of the 
previous iterations. The highest gravity coefficient correlations were now the 
S 54 S 74 and C 54 C 74 pairs (-0.998 and -0.996 respectively) with the two initial 
osculating orbital eccentricities also highly correlated (-0.997). The 66 , 86 and 
55, 75 correlations in this case were smaller in magnitude than 0.850 except for 
the S 55 S 75 correlation which was -0.926. 

Each parameter estimation run also determined the spread of the 
parameter correlations. For each different measurement type and orbital 
geometry investigated, this distribution was divided by the total number of 
correlations. Figure 7.5-1 plots this normalized distribution for the parameter 
correlations greater than 0.50. This graph shows how the single orbiter, earth- 
based sensing scheme is characterized by very high parameter correlations. If 
a subsatellite and laser ranging and Doppler instrumentation are added to this 
scheme, the correlations are driven down significantly and lunar far-side 
motion is observed. If these mission modifications are not feasible, then 
augmenting the single orbiter earth-based sensing scheme with long baseline 
interferometer measurements will also reduce the parameter correlations, 
although several of the highest correlations still remain in both cases. 

According to the graph, a better mission modification would be to 
include an elliptical "viewing" satellite to observe the near-circular low 
altitude polar "sensing" satellite. The chart suggests that if the viewing 
satellite is placed in an orbital plane orthogonal to the "sensing" satellite's 
orbital plane, the lowest parameter correlations are achieved. If both 
spacecraft, however, are launched together and then separated after Lunar 
Orbit Insertion (LOI), as was planned for the Lunar Observer mission [39], 
then coplanar "viewing" and "sensing" satellites will still provide one of the 
lowest sets of parameter correlations. 

These attempts to break the parameter correlations demonstrate that 
although some of the parameter correlations can be reduced, using a near¬ 
circular polar 200 km altitude satellite as the "sensing" vehicle to estimate an 
8 X 8 spherical harmonic expansion to model the lunar gravity field 
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consistently results in extremely high gravitational parameter correlations 
between the C 55 C 75 , £> 55575 , S 66 S 86 / and C 66 C 86 gravitational parameters. 
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Figure 7.5-1: Parameter Correlation Distributions for Several Observation Methods 

Because the highest parameter correlations are always for the n- 2 , n-2 
and n, n-2 parameters and the n-3, n-3 and n-1, n-3 parameters, the spherical 
harmonic expansion geometry should be investigated. Both sets of 
correlations involve sectorial terms and a tesseral counterpart with the same 
number of longitudinal slices. The n=m =6 sectorial and n= 8 , m =6 tesseral 
terms zero lines are illustrated by the globes in Figure 7.5-2. The geometrical 
relationship of these highly correlated gravitational parameters suggests that 
they may result from estimating the lunar gravity field solely from 
observations of a polar lunar satellite, since the satellite's ground track 
repeatedly traverses the sectorial slices of the moon. 

This analysis suggests that the best way to break the high gravity 
coefficient correlations would be to use multiple inclination "sensing" 
satellites. In his lunar harmonic gravity analysis, Alfred Ferrari noted that 
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"...[W]ell conditioned estimates of the gravity harmonics are only achieved 
when data from many different inclinations are used, ..."[19]. Multiple 
inclination observations could be achieved by using new polar satellite data 
along with existing Apollo-era near equatorial satellite data. A better scheme 
is to place near-circular, low altitude "sensing" satellites in 90° and 45° 
inclination orbits and observe their lunar far-side motion with a single 
elliptical "viewing" satellite. This data could also be combined with existing 
Apollo-era Doppler tracking data when estimating the lunar gravity field. 
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Conclusions 


8.1 Summary of Results 

This thesis investigated the ability of spherical harmonic expansion 
estimates of the lunar gravity field to predict low altitude lunar orbits 
globally. These lunar spherical harmonic expansions were estimated from 
simulated observations of a near-circular polar lunar satellite. Several 
different observation geometries and measurement methods were 
investigated; two were used to estimate the lunar gravity field. 

Since most of the methods used to derive the lunar gravity field 
employ Doppler observations, a test was performed to determine the impact 
sensing geometry had on gravity field estimation. Since the Doppler 
observations measure velocity along the line-of-sight, the geometrical 
orientation between the Doppler line-of-sight and the vector between the 
mascon and orbiting body will affect the observation's ability to detect the 
mascon's presence. For the first scenario in this test, a pair of mascons were 
placed on the lunar limbs and an eighth degree and order spherical harmonic 
expansion was estimated. For the second scenario, two mascons were placed 
on the lunar near- and far-side faces and the process was repeated. 
Unfortunately, the test case was unable to conclusively establish the impact of 
viewing geometry upon the ability to detect and model local mascon 
disturbances. 

A lunar gravity "truth" model was developed for this thesis which 
combined a fifth degree and order spherical harmonic expansion and nine 
major mass anomalies or mascons. The seven most significant lunar maria 
mass anomalies, as identified by Wong et. al. [47], served as the lunar near- 
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side mascons. The two significant far-side mascons were positioned to be 
difficult to sense from earth-based observations and to balance the lunar 
center of mass. This truth model, since it contained mascon disturbances, was 
intended to be a difficult gravity field to model with a pure spherical 
harmonic expansion. 

Using the lunar gravitational "truth" model, the near-circular polar 
lunar satellite gravity field "sensing" orbit was numerically integrated for 
fourteen days, half of a lunar period, to provide complete lunar surface 
coverage as the moon rotated under the orbital plane. Since this satellite does 
not actually measure the lunar gravity field, different techniques were used to 
observe this satellite's motion and estimate the gravity field. Estimates of the 
coefficients in spherical harmonic expansions sought to match the 
observations of the "sensing" satellite's motion. Based on the observation 
residuals between the initial parameter guesses' theoretical observations and 
the "truth" model observations, the initial conditions and harmonic 
coefficients were adjusted and the process was iterated until the best fit to the 
observations was obtained. Although different degree and order spherical 
harmonic expansion fits were attempted, the iterative process failed to 
converge for those above eighth degree and order, probably because of the 
difference between the "truth" (spherical harmonics plus mascons) and fit 
(spherical harmonic expansion only) models. 

Estimated lunar gravity fields were obtained for two different sensing 
schemes. The first scheme recreated the gravitatonal sensing method used 
during the Apollo era with mostly near-equatorial satellites. Earth-based 
Doppler observations of the near-circular polar lunar satellite were used to 
estimate the lunar gravity field. The second sensing scheme employed a 
second lunar satellite in an elliptical orbit, viewing the first satellite's motion. 
Earth- and satellite-based Doppler observables simulated the coherent bent 
pipe link between an earth tracking station, the elliptical "viewing" satellite, 
and the circular "sensing" satellite proposed for the Lunar Observer mission 
[16,39,40]. 

For the cases in which the estimation process converged upon a fit to 
the observations, the fit model was analyzed to determine how well it 
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modeled the "true" lunar gravity field. The model's fit to the observations as 
well as its ability to model global radial accelerations were investigated. 

The observation techniques which employed Doppler observations of 
the lunar far-side motion provided a much better fit to the "true" gravity 
field. Both estimation techniques employed an eighth degree and order 
spherical harmonic expansion. In neither case, however, did the eighth 
degree and order spherical harmonic expansion closely model the "true" 
lunar gravity field. 

Parameter covariance information determined in the estimation fitting 
process was used to predict satellite state uncertainties ahead one lunar orbit. 
A low altitude, low lunar inclination orbit was selected for this analysis. 
Significantly lower state covariances were predicted for the estimated lunar 
gravity field based on the dual orbiter sensing scheme. Observations of the 
lunar polar satellite during far-side passes significantly reduced the predicted 
uncertainties for the estimated fit model. In the best case, the eighth order 
spherical harmonic expansion predicted uncertainties of close to three 
quarters of a kilometer in position and one meter per second in velocity for 
the orbit, hardly acceptable for future lunar missions. 

Using the lunar gravity "truth" model, the satellite's "true" position 
and velocity were numerically integrated ahead one orbit. The true spacecraft 
state was then compared to the predicted state to reveal the state errors. These 
errors were then compared to the covariance uncertainties. Once again, the 
estimated fit model which included lunar far-side observations did a 
significantly better job of predicting the correct satellite state, although it still 
provided errors as large as 2.8 kilometers in position and 2.3 meters per 
second in velocity. In neither case did the spherical harmonic expansion fit 
model's covariance uncertainties come close to predicting the correct state 
errors, further evidence of the eighth degree and order expansion's inability 
to successfully predict lunar orbits for future lunar missions. 

The estimated fit models were then used to plan a lunar landing 
deorbit maneuver. Starting on the near side of the moon in a circular low 
lunar inclination parking orbit, the satellite's orbit was numerically integrated 
to the far side of the moon where a deorbit maneuver placed it in an elliptical 
transfer orbit. When the spacecraft reached a specified lunar altitude, as 
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detected by an on-board radar altimeter, it began its powered descent to the 
lunar surface. After planning the deorbit maneuver to reach a specific 
location for powered descent initiation (PDI), the deorbit maneuver was 
executed in the "true" lunar gravity field. 

The lunar deorbit maneuver planned with the single orbiter earth- 
based Doppler sensing scheme estimated gravity field reached PDI altitude 
eight degrees in selenographic longitude and three quarters of a degree in 
latitude away from its target position, a surface error of over two hundred and 
forty kilometers. The lunar deorbit maneuver planned for the dual orbiter 
gravity field was only one and three quarters of a degree of longitude and one 
fifth of a degree of latitude off target when it reached PDI altitude. This 
scenario's lunar surface error was fifty-six kilometers, an error much easier to 
recover from during powered descent. Both cases, however, indicate that 
further navigation aides will be required during lunar operations, such as 
radio beacons at a lunar base or a global lunar navigation system. 

Extremely high parameter correlations were encountered when 
estimating the gravity coefficients. Three different measurement schemes 
were investigated to see if they could reduce these parameter correlations. In 
the first case, the elliptical "viewing" satellite was placed in an orbital plane 
orthogonal to the circular "sensing" satellite's orbital plane. In the second 
case, earth-based interferometer measurements were added to the single 
orbiter, earth-based Doppler sensing scheme. Both east-west and north-south 
look angles were simulated in this search for a way to break the high 
parameter correlations. The third case simulated the Goddard Space Flight 
Center's proposed co-orbital laser ranging and Doppler sensing scheme [2]. 

The first method was able to reduce some of the high parameter 
correlations from the dual orbiter co-planar sensing scheme, but several of 
the highest parameter correlations remained. The second case significantly 
reduced the high parameter correlations from the single satellite earth-based 
Doppler-only sensing scheme, although once again the same highest 
parameter correlations still remained. The third case produced some 
different high parameter correlations. This method reduced the parameter 
correlations significantly from the single orbiter, earth-based Doppler sensing 
scheme. Although this case would reduce the parameter correlations from 
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the single orbiter case, it does not reduce the correlations as much as the dual 
orbiter coplanar and orthogonal cases. 

Since none of these new observation techniques eliminated the highest 
parameter correlations, the spherical harmonic expansion geometry was 
investigated to determine why the same parameters were so highly 
correlated. It seems that the extremely high correlations result from using a 
single polar "sensing" orbit for all of the gravity field estimates. Polar ground 
tracks made it difficult to separate the effects of some of the sectorial and 
related tesseral terms. "Sensing" orbits at different orbital inclinations should 
break the extremely high parameter correlations and aid the lunar gravity 
field estimation process. 


8.2 Recommendations for Future Research 

The research for this thesis revealed several topics for further 
investigation. First, a more efficient gravity field model based on a surface 
layer representation should be investigated. A surface layer representation 
would require roughly one third of the parameters than a spherical harmonic 
expansion to model the lunar gravity field's high frequency behavior. The 
surface layer representation should be able to constrain the total lunar mass 
and lunar center of mass. Estimated fits with this gravity field model, if it is 
programmed into the Planetary Ephemeris Program, can be compared to 
those obtained with a spherical harmonic expansion gravity field model. 

Additionally, this thesis only compared the estimated fit models for 
two observation techniques. Because of its observability into the lunar far- 
side motion, the dual orbiter scheme naturally provided a better fit to the 
lunar gravity "truth" model. Estimated fit models should be developed for 
the three different observation methods investigated to break the high 
parameter correlations. Different Doppler observation methods, besides the 
bent pipe method, can be investigated, including the satellite bounce and 
satellite beacon methods proposed by JPL [40]. Straightforward software 
changes in PEP should allow simulation of these observables, and further 
modifications would allow the incorporation of lunar navigation aids (either 
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lunar surface sites or navigation satellites). Additionally, two "sensing" 
satellites in different inclinations should be studied to verify their ability to 
break the high parameter correlations. Higher degree and order simulated fits 
should also be attempted to determine the navigation improvements that can 
be achieved with higher resolution spherical harmonic expansions. 

For further studies, a less stringent lunar gravity "truth" model should 
be used. The difficult to observe mascon on the lunar far side was too large, 
in both surface area and total strength. The lunar center of mass should still 
be constrained, but the far-side mascons should be roughly the same size and 
strength as their near-side counterparts. 

If lunar rotation moment of inertia partial derivatives in the Planetary 
Ephemeris Program are changed to equivalent second harmonic partial 
derivatives, PEP can include lunar laser corner reflector observations 
simultaneously with satellite observations in the estimation of the lunar 
gravity field. The lunar laser observations would provide very accurate 
determinations of the lower degree harmonics (second and even third and 
fourth degree [13]), while the lunar satellite observations would aid the 
estimation of higher degree harmonics. 

Finally, the most appropriate method to estimate the lunar gravity 
field with real observations is to use maximum likelihood system 
identification. This method runs a Kalman filter on the satellite motion with 
noise in the dynamics due to unmodeled forces, and applies a maximum 
likelihood estimator to gravity, initial condition, and other parameters [31 
Chapter 10]. 
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Evaluation of Legendre Polynomials 

and 

Normalized Legendre Functions 


A.1 Legendre Polynomials 

The Legendre polynomials and their first and second derivatives are 
required for the numerical integration of the equations for satellite motion 
and for the partial derivatives of satellite motion. The Legendre polynomials 
are used to determine the zonal harmonic gravity effects and are also used in 
the recursive formulas for the Legendre functions. Unnormalized zonal 
harmonic gravity coefficients are input to PEP since its algorithms use the 
unnormalized form of the Legendre polynomials. These algorithms were left 
unchanged, because the normalization factor in the equation 

P„(z) = V2n + l P„(z ) (A. 1-1) 

does not grow too rapidly with zonal harmonic degree n, and the normal 
equation's automatic scaling feature works around any numerical problems 
(Section 5.1.6). A switch to normalized zonal harmonic coefficients C n p from 
J n would also involve a sign change in several PEP subroutines. 

The recursive evaluation of P n (z), Pn(z), Pn(z) in subroutine LEGNDR 
(and the new subroutine LEGNDS) uses the relations [1, 8] 

nP„ (z) = (2« -1) z P„_, (z) - (n - 1)P„_2 (z) (A.l-2) 

P'„ (z) = P;_ 2 (z) + (2 n -1 )P„_, (z) (A.l-3) 

P''(z) = P" 2 (z) + (2 n- 1)P;_, (z) (A. 1-4) 
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with the starting values 

P 0 (z) = 1, P, (z) = z , P 2 (z) = %(3z 2 -1), P 3 (z) = K(5z 3 -3z) (A.l-5) 

P'(z) = 0, P,'(z) = l, Pj(z) = 3z, P^(z) = K(15z 2 -3) (A.l-6) 

P"(z) = 0, P"(z) = 0, Pj(z) = 3, P"(z) = 15z (A.1-7) 

The evaluation of the above functions starts with n = 2, since the 
spherical harmonic expansions in PEP start at the second degree (assuming 
the center of mass of the central body coincides with the origin of 
coordinates). 


A.2 Normalized Legendre Functions 


The Legendre functions and their first and second derivatives are 
required for the numerical integration of the equations for satellite motion 
and for the partial derivatives of satellite motion. Legendre functions are 
used to determine the tesseral harmonic gravity effects. The original version 
of PEP converted normalized tesseral harmonic gravity coefficients (C nm/ 
S nm ) from the input stream into unnormalized coefficients (C nm / S nm ) since 
its internal algorithms computed the unnormalized versions of the Legendre 
functions. The option was added to PEP (incorporated into the SAO version) 
to use normalized Legendre functions in the numerical integration process. 
This option was desirable since the normalization factors vary widely for high 
degrees n, especially as m approaches n in the equation 


P„Jz) = 


2(2n + l)(n-m)! 
V (n + m)\ 




m = l,...,n (A.2-1) 


Define the precalculated coefficients for m = l,...,n and n > 1 


a 


nm 


(2w + l)(n-m) 
(2n-l)(n + m) 


(A.2-2) 
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l(2n)(2n + l) 

m = 1 


V <w + D ' 


(A.2-3) 

l(2n + l)(n + m-l) 

1< m< n 


V (2n-l)(n + m) 



= V(2n)(2n + l)(n + l) , 

m - 1 

(A.2-4) 

= ^(n + m)(n-m + 1), 

1 < m < n 



and the following expression from the recursive formulas for the 
unnormalized Legendre functions 

q = Vl-z 2 (A.2-5) 

The recursive evaluation of P nm (z), P^z), P; m (z) is performed in the new 
subroutine LEGNDS using the relations below [1, 8]. The recursive formulas 
for the normalized Legendre functions are 

P n (z) = V3^, P 21 (z) = Vl5 q z , P 22 (z) = ^Vl5 q 2 (A.2-6) 

P n l^ = a n\ 2 ^-1,1 ^ *7 Pn-l( Z ) ' (A.2-7) 

+ ™ = 2 .< A - 2 ' 8) 

Pnn & = b m q (A.2-9) 

The recursive formulas for the first derivative of the Legendre functions with 
respect to the argument z (using unnormalized Legendre polynomials) are 

P'(z) = —P 2 ',(z) =—(l-2z 2 ), P 2 ' 2 (z) = -zVl5 (A.2-10) 

q q 

P' t (z) = f (z P nl (z) - d„, q P n (z)) (A.2-11) 

P; m (z) = £(mzP n Jz)-d„ m c,P n , m ^z)), m = 2,...,n (A.2-12) 

Lastly, the recursive formulas for the second derivative of the Legendre 
functions with respect to the argument z (using unnormalized Legendre 
polynomials) are 
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P"( z) = ~, P 2 ';(z) = ^( 2z 2 -3), P£(z) 


= -Vl5 


p:, m=\ 


( (J -(- 2^)_ _ ^ 

ii±fiP nl (2) + 2P' 1 (2) 

-rf„,(? P ” (2) + ‘7^ <z) ) 


(A.2-13) 


(A.2-14) 



f f (1 + 2 2 )P nm (z) + m 2 P' nm ( 2 ) ^ 
(f <Z) + <7 ?n,*-l < Z >) V 


m = 2,...,n 

(A.2-15) 


The recursive evaluation of the above functions starts with n = 2 for 
the gravity harmonic force evaluation, since the spherical harmonic 
expansions in PEP start at the second degree (assuming the center of mass of 
the central body coincides with the origin of coordinates). 
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Appendix B 


List of Constants 


Primary 

Constant 

Description 

Value 

Units 

[20] 

M s 

Mass of the sun 

i 


-Jgm. 

Gaussian gravitational constant 

0.01720209895 

AU3/2/day 

CT 

Coordinate Time 

A.l + 32.15 s 
TAI+32.184 5 

day 

A.l, TAI 

Atomic Time Seconds 

9,192,631,770 

oscillations 
of cesium 


The A.l and TAI atomic times are kept as the average of a number of cesium atomic 
clocks at the national time services (particularly at the U.S. Naval Observatory). UTC time 
runs at the A.l and TAI rates, and presently differs from TAI by an integral number of seconds. 
Every six or twelve months there is an increase of one second (a UTC leap second) in TAI-UTC, 
to keep UTC within one second of UT1 time, defined by a formula in terms of sidereal time [23]. 
Given a UTC observation time, CT is computed for interpolating from ephemeris files, and UT1 
time and earth wobble are computed from tables published by the U.S. Naval Observatory and 
the Bureau International de l'Heure (see Appendix C.2). 
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Auxiliary 

Constant 

[23] 

Description 

Value 

Units 

AULTSC 

Astronomical Unit 

in Light Seconds 

499.004782 

seconds/AU 

c 

Speed of light 

299,792.458 

km/sec 

M s 

M c 

sun/(earth + moon) mass ratio 

328,900.1 


M c 

M, 

(earth + moon)/moon mass ratio 

82.301 


GM, 

Lunar gravitational constant 
(PEP uses AU 3 /day 2 ) 

4,902.79375 

km 3 /sec 2 

Pi 

Lunar mean equatorial radius 

1,738 

km 

P / 

Lunar period 

27.322 

days 

gal 

Unit of Acceleration 

1.0 

cm/ sec 2 
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Inertial Coordinate Transformations 


C.1 Transformation Between Moon Fixed and Inertial Frames 

The transformation between moon-fixed (selenographic) coordinates 
and the inertial integration frame coordinates (the mean equinox and equator 
of 1950.0) is required to evalute gravity harmonic and mascon accelerations in 
the numerical integration of a lunar satellite's motion. Auxiliary software 
also uses this transformation to analyze integration output in the 
selenographic frame and to prepare numerical integration input. In Equation 
(3.2-2) for the transformation between coordinates fixed in the moon and the 
integration frame coordinates, the rotation matrix R can be expressed as [6] 

R = UVP (C.l-1) 


where 


P = Earth precession matrix transforming between integration 
frame coordinates and coordinates referred to the mean 
equinox and equator of date (50 arcseconds per year). 

V = Transformation between coordinates referred to the mean 
equinox and equator of date and coordinates referred to the 
mean equinox and ecliptic of date (23.4° rotation). 

U = Transformation between coordinates referred to the mean 
equinox and ecliptic of date and coordinates fixed in the 
moon along the nominal principal moment of inertia axes 
(360° per 27.2 days). 
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PEP uses a Taylor series expansion to calculate the earth precession 
matrix (P) for speed during a numerical integration [7]. If e 0 is the standard 
expression for the obliquity of the ecliptic [20, 23], the transformation matrix 
between mean equatorial and ecliptic coordinates of date is 

0 0 

cos £q sine 0 (C.l-2) 

- sin £q cose 0 

Cassini's laws plus the physical libration of the moon determine the 
transformation between the mean equinox and ecliptic of date and 
selenographic coordinates (U). The following formulas are used within PEP 
to calculate this transformation matrix. 

Letting 

M = Mean longitude of the moon measured in the ecliptic from 
the mean equinox of date to the mean ascending node of 
the lunar orbit and then along the orbit (27.2° day period). 

£2 = Longitude of the mean ascending node of the lunar orbit 

on the ecliptic measured from the mean equinox of date 
(18.6 year period). 

I = Inclination of the lunar equator to the ecliptic (1.53889°). 

a, p, x = Physical librations in node, inclination, and longitude. 

Standard polynomial expressions are used for the angles M and £2, as well as 
for the angles £, V, F, and D which are from Brown’s lunar theory [23]. PEP 
uses the following trigonometric series for the physical librations [28] 

x = - 12.9" sin i - 0.3" sin 21 + 65.2" sin V + 9.7" sin (2f - 21) 

+ 1.4" sin (2F - 2D) +2.5" sin (D - £) - 0.6" sin (2D - 21 + £') 

- 7.3" sin (2D - 2£) - 3.0" sin (2D - t) - 0.4" sin 2D 

+ 7.6" sin £2 (C.l-3) 

p = -106" cos t + 35" cos (2F - £) - 11” cos 2F 

- 3" cos (2F - 2D) - 2" cos (2D - £) (C.l-4) 
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I(x - a) = 108" sin t - 35" sin (2F - l) + 11" sin 2F 

+ 3" sin (2F - 2D) + 2" sin (2D - () (C.l-5) 

where I is measured in radians. 

In terms of these angles, the Euler angles for the moon rotation are 

y/ = £2 + <7 

d = l+p (C.l-6) 

<p = 180 o +(M-£2)+(t-<7) 

The transformation matrix between selenographic coordinates fixed in the 
moon and those referred to the mean equinox and ecliptic of date is then [ 6 ] 

Un = -siny cos 0 sin (j) - cosy cos<|> 

U ]2 = cos y cos 0 sin <(> + sin y cos <|> 

U ]3 = - sin 0 sin <(> 

U 21 = - sin y cos 0 cos <J) - cos y sin <|> 

U 22 = cos y cos 0 cos - sin y sin <|) 

U 23 = - sin 0 cos <t> 

U 31 = - sin y sin 0 

U 32 = cos y sin 0 

U 33 = cos 0 (C.l-7) 

To check these transformations, the earth and sun's selenographic 
longitude and latitude were printed out for certain dates. These positions 
were obtained by transforming the vectors pointing from the moon to these 
objects with the transformation matrix Equation (C.l-1). This check verified 
that PEP's values agreed with the published values in the Astronomical 
Almanac, to the number of places published (0.001° for the earth, 0.01° for the 
sun). 

Since this transformation's moon-fixed axes are not exactly along the 
moon's principal moments of inertia axes, the C 21 , S 21 , and S 22 spherical 
harmonic coefficients cannot be assumed to be zero. If the rotation and orbit 
of the moon were estimated by fitting to lunar laser corner reflector data, 
these second degree harmonic coefficients could be set to zero. 
Simultaneously processing lunar laser data with lunar orbiter data could 
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provide the best estimate of the remaining second degree harmonic 
coefficients, as well as third and perhaps fourth degree coefficients. 

One auxiliary program (selenang), converts satellite osculating orbital 
angles in the selenographic coordinate frame to angles in the integration 
(1950.0) frame. These angles are input to PEP to specify a desired 
selenographic orbital orientation for numerical integrations. Within PEP, 
these initial osculating elliptic orbital elements are converted to Cartesian 
position and velocity initial conditions for the numerical integrations. 

Another auxiliary program (selenelm) converts the numerical 
integration output to selenographic position and velocity and to 
selenographic osculating elliptic orbital elements as functions of time. The 
selenographic osculating elliptic orbital elements, the altitude above the lunar 
surface, and the subsatellite selenographic longitude and latitude can then be 
plotted over time. 


C.2 Transformation Between Earth Fixed and Inertial Frames 

Since earth-fixed observing sites need to be transformed to the inertial 
frame to process observations, the transformation matrix between earth-fixed 
coordinates and the inertial integration frame referred to the mean equinox 
and equator of 1950.0 is also required. This rotation transformation matrix R 
can be expressed as [8] 

R = WSNP (C.2-1) 


where 

P = Earth precession matrix transforming between integration 
frame coordinates and coordinates referred to the mean 
equinox and equator of date (50 arcseconds per year). 

N = Earth nutation matrix transforming between coordinates 
referred to the mean equinox and equator of date and 
coordinates referred to the true equinox and equator of date 
(20 arcseconds). 
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S = Earth sidereal time matrix transforming between 
coordinates referred to the true equinox and equator of date 
and coordinates with z axis along the earth pole of rotation 
and x axis in the meridian of Greenwich through the 
rotation pole (360° per 23 h 56 m 4.09 s ). 

W = Earth wobble matrix transforming between coordinates 
with z axis along the earth pole of rotation and x axis in the 
meridian of Greenwich and coordinates with z axis along 
the conventional international pole and x axis in the 
meridian of Greenwich through the conventional 
international pole (0.3 arcseconds). 

Since PEP uses the IAU value of the precession constant [23] in the 
trigonometric angles for evaluating the precession matrix, P, numerical 
integration results in the 1950.0 reference frame are transformations of 
integration results in the IAU J2000.0 reference frame [23]. 

Since observations in PEP are a function of UTC time and there is a 
mathematical formula relating sidereal time to UT1 time [23], the 
relationship between UTC and UT1 time is needed to determine the sidereal 
time transformation matrix, S. PEP determines this relationship from values 
published by the U.S. Naval Observatory based on photographic zenith tube 
and other observations. 

UT1 - UTC = (A.l - UTC) - (A.1 - UT1) (C.2-2) 

The U.S. Naval Observatory and the Bureau International de l'Heure 
also publish wobble coordinates based on the above mentioned observations 
which are used to compute the earth wobble transformation matrix, W. Exact 
values of these quantities were not required for this thesis' simulations. 
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Appendix D 


Lunar Gravitational "Truth" Model 


Table D-l: 5X5 Spherical Harmonic Coefficient Portion of Truth Model [121 


Harmonic 

xlO- 6 

Harmonic 

x 10" 6 

h 

202.431 



C 2 i 

-0.07 

S 21 

-0.00 

C 2 2 

34.49 

S 22 

0.03 

h 

8.8897 



C31 

21.96 

S31 

6.63 

C 32 

14.14 

S 32 

4.76 

C 33 

15.87 

S 33 

-2.45 

U 

-11.73 



C 41 

-4.82 

S 41 

1.91 

C 42 

-8.13 

S 42 

-6.76 

C 43 

0.48 

S 43 

-14.43 

C 44 

-3.50 

S 44 

-0.55 

J 5 

2.388 



C 51 

-9.66 

S 51 

-1.53 

C 52 

3.71 

S 52 

-2.35 

C 53 

-0.39 

S 53 

4.91 

C 54 

0.56 

S 54 

-6.58 

C 55 

-6.69 

S 55 

11.60 
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Table D-2: Lunar Near-Side Mascons at 1638 km Radius [47] 


Mass 

East 

North 

Mass 

East 

North 

(10^ Lunar 

Longitude 

Latitude 

(10“^ Lunar 

Longitude 

Latitude 

Mass) 

(degrees) 

(degrees) 

Mass) 

(degrees) 

(degrees) 

Sea of Serenity 

| Sea of Rains 

0.101 

10.5 

34.0 

-0.597 

342.0 

48.0 

2.437 

17.5 

34.0 

0.467 

328.5 

40.5 

2.453 

10.5 

28.0 

5.329 

337.5 

40.5 

3.316 

17.5 

28.0 

5.218 

346.5 

40.5 

3.774 

24.5 

28.0 

1.252 

328.5 

34.0 

0.371 

8.25 

22.5 

1.628 

335.5 

34.0 

2.628 

13.75 

22.5 

3.586 

342.5 

34.0 

3.628 

19.25 

22.5 

1.752 

349.5 

34.0 

1.681 

24.75 

22.5 

-0.421 

328.5 

28.0 

0.290 

12.5 

17.5 

2.336 

335.5 

28.0 

0.624 

17.5 

17.5 

2.654 

342.5 

28.0 

0.229 

22.5 

17.5 

Seething Ba 


| Sea of Crises | 

■ 


12.5 

0.230 

52.25 

22.5 



12.5 

1.520 

57.75 

22.5 

0.531 

| 

7.5 

0.992 

52.5 

17.5 

1.522 

352.5 

7.5 

2.837 

57.75 

17.5 

Sea of Moisture 

1.836 

62.5 

17.5 

3.091 

318.75 

-22.5 [ 

1.183 

57.5 

17.5 

2.726 

325.25 

-22.5 

0.688 

62.5 

12.5 

0.164 

321.5 

-28.0 

| Sea of Nectar j 

Smyth's Sea 

■EEEfH 

27.5 

-12.5 

0.622 

87.5 

2.5 

1.394 

32.5 

-12.5 

-0.172 

82.5 

-2.5 

1.199 

37.5 

- 12.5 

0.873 

87.5 

-2.5 


27.5 

-17.5 

0.739 

92.5 

-2.5 | 


32.5 

-17.5 

0.352 

87.5 

-7.5 


37.5 

-17.5 

0.691 

92.5 

-7.5 
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Table D-3: Lunar Far-Side Mascons at 1638 km Radius 


Mass 

East 

North 

Mass 

East 

North 

(10 -6 Lunar 

Longitude 

Latitude 

(10- 6 Lunar 

Longitude 

Latitude 

Mass) 

(degrees) 

(degrees) 

Mass) 

(degrees) 

(degrees) 

Difficult to Observe Mascon 

0.829 


8.0 

1.836 


-8.0 

1.504 

180.0 

8.0 

2.367 


-8.0 

1.254 

177.5 

8.0 

1.726 

182.5 

-8.0 

0.615 

183.75 

5.0 

1.357 


-5.0 

1.391 

181.25 

5.0 

2.929 

KKraEZH 

-5.0 

1.613 

178.75 

5.0 

2.265 

181.25 

-5.0 

0.825 

176.25 

5.0 

0.935 

183.75 

-5.0 

2.125 

185.0 

0.0 

1.323 

175.0 

0.0 

4.737 

182.5 

0.0 

2.682 

177.5 

0.0 

7.712 

180.0 

0.0 

Balancing Mascon 

1.602 

WUdsM| 

-50.5 

4.151 


-48.0 

2.534 


-50.5 

7.212 


-48.0 

5.128 

213.25 

-45.5 

2.817 


-48.0 

3.713 

215.75 

-45.5 

Final Point Masses 

1.11385 

217.6253 

Bll 


214.93 

-49.06 

2.1xl0 -11 

203.77 

-43.54 
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Appendix E 


Tables of Spherical 
Harmonic Coefficients 


On the following pages are the tables containing the spherical 
harmonic coefficients for the models estimated in this thesis. Also included 
are the spherical harmonic coefficients from degree six through ten from Alex 
Konopliv's 50 X 50 spherical harmonic model estimated at the Jet Propulsion 
Laboratory. These coefficients were used in the fit of a 10 X 10 model to 
observations generated by a 10 X 10 truth model. The tesseral coefficients are 
normalized and the zonal coefficients are unnormalized to follow the 
convention used in PEP-D which was modified from the SAO's version to 
use normalized Legendre polynomials, but not normalized Legendre 
functions. 



’fiLn? 


'•hr* 
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Table E-l: Limb Mascon Fit Model Coefficients 



Note: All of the C nm and Snm terms are Normalized. The J n terms are Unnormalized. 
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Table E-2: Face Mascon Fit Model Coefficients 


Harmonic 

xlO' 6 

Harmonic 

xlO' 6 

h 

212.9685 


-18.2734 

h 

-17.2629 

■ mm 

56.7533 

J6 

6.8783 

J7 

-33.1132 

J8 

1.9000 



c 2 i 

10.1197 

S 2 1 

■■jjlglg 

C22 

35.1938 

S22 


C31 

20.1172 


8.9737 

C32 

0.1584 


5.8036 

Q3 

19.8874 


-5.5137 

Qi 

-25.7774 


3.4003 

Q2 

-2.1818 


-11.9870 

C43 

14.6293 


-18.4118 

C44 

-4.5549 

S44 

3.9231 

C51 

-7.5438 

S51 

-0.8299 

C52 

24.1943 

S52 

-7.5913 

C53 

-8.1387 

S53 

12.9983 

C54 

-11.5970 

S54 

-0.2905 

C55 

-4.6557 

S55 

7.4909 

C61 

17.6456 

S6i 

-3.3481 

Ol 

-0.9812 

S62 

1.0038 

03 

-16.1761 

S63 

7.6098 

Q>4 

4.9624 

S64 

-9.3866 

05 

6.4501 

S65 

-4.1229 

06 

-0.9473 

S66 

3.4506 

C71 

1.5589 

S71 

-3.8029 

C72 

-9.6129 

S72 

4.0109 

C73 

2.9190 

S73 

-4.8617 

C74 

10.2397 

S74 

-4.9557 

C75 

-3.1597 

S75 

7.1339 

C76 

-0.5898 

S76 

1.1869 

C77 

0.1187 

S77 

0.0073 

Q?i 

-6.5178 

S81 

-0.4581 

C82 

-2.6772 

S82 

1.4211 

03 

3.3579 

S83 

-4.3328 

04 

-0.5695 

S84 

4.2547 

05 

-5.7318 

S85 

2.0067 

06 

3.6244 

S86 

-5.2175 

07 

-0.4581 

S87 

-0.4756 

08 

1.5694 

S88 

0.0024 



Note: All of the C n m and Snm terms are Normalized. The J n terms are Unnormalized. 
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Table E-3: 8X8 Single Orbiter Earth-based Doppler Fit Model Coefficients 



Note: All of the C n m and S nm terms are Normalized. The J n terms are Unnormalized. 
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Appendix E: Tables of Svherical Harmonic Coe 


Table E-4: Additional Coefficients for the 10 X 10 Spherical Harmonic Expansion [27] 



Note: All of the C n m and S nm terms are Normalized. The J n terms are Unnormalized. 
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Table E-5: 8X8 Dual Orbiter Bent Pipe Doppler Fit Model Coefficients 



Note: All of the C nm and S nm terms are Normalized. The J n terms are Unnormalized. 
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